A new method for measuring azimuthal distributions in nucleus-nucleus collisions 



O 

o 

(N 

U 

Or 
< 



X 



Nicolas Borghini, 1, * Phuong Mai Dinh, 2 and Jean- Yves Ollitrault 2 

1 Laboratoire de Physique Theorique des Particules Elementaires 
Universite Pierre et Marie Curie, 4 place Jussieu, F-75252 Paris cedex 05 
2 Service de Physique Theorique, CEA-Saclay, F-91191 Gif-sur-Yvette cedex 

The methods currently used to measure azimuthal distributions of particles in heavy ion colli- 
sions assume that all azimuthal correlations between particles result from their correlation with the 
reaction plane. However, other correlations exist, and it is safe to neglect them only if azimuthal 
anisotropies are much larger than 1/y/N, with N the total number of particles emitted in the colli- 
sion. This condition is not satisfied at ultrarelativistic energies. We propose a new method, based on 
a cumulant expansion of multiparticle azimuthal correlations, which allows measurements of much 
smaller values of azimuthal anisotropies, down to 1/iV. It is simple to implement and can be used to 
measure both integrated and differential flow. Furthermore, this method automatically eliminates 
the major systematic errors, which are due to azimuthal asymmetries in the detector acceptance. 



I. INTRODUCTION 



In heavy ion collisions, much work is devoted to the study of the azimuthal distributions of outgoing particles, 
and in particular of distributions with respect to the reaction plane. Since these distributions reflect the interactions 
between particles, possible anisotropies, the so-called "flow," reveal information on the hot stages of the collision: 
, thermalization, pressure gradients, time evolution, etc. Q 

■ Since the orientation of the reaction plane is not known a priori, flow measurements are usually extracted from 
^T) \ two-particle azimuthal correlations. This is based on the idea that azimuthal correlations between two particles are 

• generated by the correlation of the azimuth of each particle with the reaction plane. The assumption that this is the 
only source of two-particle azimuthal correlations, or at least that other sources can be neglected, dates back to the 
early days of the flow It still underlies the analyses done at ultrarelativistic energies, both at the Brookhaven 
AGS and the CERN SPS. 

However, we have shown in recent papers ||,|| that other sources of azimuthal correlations (which we refer to as 

■ "nonflow" correlations) are of comparable magnitude and must be taken into account in the flow analysis. We have 
studied in detail the well-known correlations due to global momentum conservation Q and those due to quantum 
correlations between identical particles ||. We have also discussed other correlations due to resonance decays and 

O | final state interactions . 

Nonflow correlations scale with the total multiplicity N like like 1/N. Thus, they become large for peripheral 
collisions. It is important to take them into account, in particular, when studying the centrality dependence of the 
flow, which has been recently proposed as a sensitive probe of the phase transition to the quark-gluon plasma |^| ^j. 

Clearly, a reliable flow analysis should eliminate nonflow correlations. Correlations due to momentum conservation 
can be calculated analytically and subtracted from the measured correlations, so as to isolate the correlations due 
to flow short range correlations can be measured independently and subtracted in the same way |^|. Other 

well-identified nonflow correlations can be estimated through a Monte-Carlo simulation ||. This method was used 
by the WA93 Collaboration to estimate direct correlations from n° — > 77 decays |l0| . Alternatively, one can attempt 
to eliminate nonflow correlations directly at the experimental level: effects of momentum conservation cancel if the 
detector used in the flow analysis is symmetric with respect to midrapidity || ; short range correlations are eliminated 
if one correlates two subevents separated by a gap in rapidity. This is the method recently used by the STAR 
Collaboration at RHIC jjlj ]. In the STAR paper, the correlations between pions of the same charge are also compared 
with correlations between ir + and 7r~: correlations from p° — > 7r + 7r _ are thus found to be negligible. 

Nevertheless, nonflow correlations remain which cannot be handled so simply. Correlations due to resonance decays, 
for instance, are hard to estimate (this would require a detailed knowledge of the collision dynamics) and there is no 
known systematic way to eliminate them at the experimental level; more importantly, the production of minijets will 
contribute to azimuthal correlations in the experiments at higher energies, at RHIC and LHC. Finally, the existence 
of other sources of nonflow correlations, so far unknown, cannot be excluded. 



*Present address: Service de Physique Theorique, CP225, Universite Libre de Bruxelles, B-1050 Brussels 



1 



The purpose of this paper is to propose a new method for the flow analysis which requires no knowledge of nonflow 
correlations. The general idea is to eliminate these latter using higher order azimuthal correlations. Higher order 
correlations were previously used in |L2| to show qualitatively the collectivity of flow. The study presented in this 
paper is more quantitative : by means of a cumulant expansion, we are able to extract the value of the flow from 
multiparticle correlations. The method we propose is more reliable, and in many respects simpler than traditional 
methods In particular, detector defects, which must be considered carefully when measuring anisotropics of a few 
percent, can be corrected in a compact and elegant way. 

In Sec. H we give the principle of our method as well as orders of magnitude. We show in particular that this 
method is more sensitive: it allows measurements of azimuthal anisotropies down to values of order l/M, instead of 
1/vN with the standard analysis, where N denotes the total multiplicity of particles emitted in the collision. 

Then, we show how the method can be implemented practically. As usual, the measurement of azimuthal distri- 
butions is performed in two steps. First, one reconstructs approximately the orientation of the reaction plane from 
the directions of many emitted particles, and one estimates the statistical uncertainties associated with this recon- 
struction. In fact, this first step amounts to measuring the value of the flow, integrated over some region of phase 
space (corresponding typically to a detector). We show in Sec. [II how this measurement can be done using moments 
of the distribution of the Q-vector, which generalizes the transverse momentum transfer introduced by Danielewicz 
and Odyniec in order to estimate the azimuth of the reaction plane Q . We also discuss an improved version of the 
subevent method introduced by the same authors to estimate the accuracy of the reaction plane reconstruction. 

The second step in the flow analysis is to perform more detailed measurements of azimuthal distributions, for 
various particles, as a function of rapidity and/or transverse momentum. We refer to these detailed measurements as 
to "differential flow" . They are usually performed by measuring distributions with respect to the reconstructed reaction 
plane, and then correcting for the statistical errors in this reconstruction, which have been estimated previously. Here, 
the differential flow will be extracted directly from the correlation between the azimuths of the outgoing particles 
and the Q-vector, as explained in Sec. IV. The discussion applies so far to an ideal detector. A general way of 
implementing acceptance corrections adapted to our method is discussed in Sec. |v|. Finally, the correct procedure is 



summarized in Sec. VI. Readers already familiar with flow analysis and willing to apply our method may go directly 
to this last section. 



II. CUMULANT EXPANSION OF AZIMUTHAL CORRELATIONS 



As the standard methods o f flow analysis (9) , ou r me thod is based on a Fourier expansion of azimuthal distributions 



L3l which is defined in Sec. II A 



Then, in Sec. II B, we discuss two-particle azimuthal correlations, on which the 
standard flow analysis relies, and show that they decompose into a contribution from flow and an additional term of 
order 1/iV which corr esponds to nonflow correlations; this latter contribution limits the sensitivity of the t raditional 



method. In Sec. II C, the decomposition is generalized to multiparticle correlations. Finally, in Sec. II D, we show 



that this decomposition of multiparticle correlations allows us to obtain more sensitive measurements of flow. 



A. Fourier coefficients 



We call "flow" the azimuthal correlations between the outgoing particles and the reaction plane. These are con- 
veniently characterized in terms of the Fourier coefficients v n |13| which we now define. In most of this paper, we 
shall work with a coordinate system in which the x axis is the impact direction, and (x, z) the reaction plane, while 
(f> denotes the azimuthal angle with respect to the reaction plane. In this frame, the momentum of a particle of mass 
m is 

(Px = Pt cos <p \ 
p y = p t sin cjy , (1) 

Pz = \/p T + m 2 sinh y / 

where pr is the transverse momentum and y the rapidity. Since the orientation of the reaction plane is unknown in 
experiments, so is the azimuth tf>. Therefore p x and p y are not measured directly. 

When necessary, we shall denote by </> the azimuthal angle in the laboratory frame. Unlike <fi, (j> is a measurable 
quantity, related to <f> by <j) = <p + <pn, where 4>r is the unknown azimuthal angle of the reaction plane in the laboratory 
system. 

With these definitions, v n can be expressed as a function of the one-particle momentum distribution /(p) = dA^/d 3 p 
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/ e i ^/(p)d 3 p 

v n (V) ee (e**) = J -2- , (2) 

/ /(P)d 3 p 
Jv 

where the brackets denote an average value over many events, and V represents a phase space window in the (pt,u) 
plane where flow is measured, typically corresponding to a detector. Since the particle source is symmetric with 
respect to the reaction plane for spherical nuclei, (sin n<j)) vanishes and v n is real. 

The purpose of the flow analysis is to extract v n from the data. Only the first two coefficients v\ and vi have been 
published. They are usually called directed and elliptic flow, respectively. There are so far very few measurements of 
higher order coefficients. The E877 experiment at the Brookhaven AGS reported values compatible with zero for V3 
and V 4 [14|. Nonvanishing values of higher harmonics, up to vq, were reported from preliminary analyses at the CERN 
SPS |15|h|. However, the latter results are likely to be strongly biased by quantum two-particle correlations ||. At 
the energies of the CERN SPS, V\ and v 2 are of the order of a few percent |l7j |, close to the limit of detectability with 
the standard methods, hence the need for a new, more sensitive method. 



B. Two-particle correlations 

Since the actual orientation of the reaction plane is not known experimentally, one can only measure relative 
azimuthal angles between outgoing particles. The standard flow analysis relies on the measurement of two-particle 
azimuthal correlations, which involve the two-particle distribution /(pi,P2) = dA/d 3 pid 3 p2: 

/ e m ^ 1 -^V( Pl ,p 2 )d 3 p 1 d 3 p 2 

/ /(Pi,P2)d J pid a p 2 

JT> 1 XT<2 



The standard analysis neglects nonflow correlations. Under that assumption, the two-particle momentum distribution 
factorizes: 

/(Pi,p 2 ) = /(pi)/(p 2 ). (4) 

Then, Eqs. (|) and (§) give 

/ e i»(*x-*)\ = Vn (p 1 )v n (p 2 ). (5) 

This equation means that the only azimuthal correlation between two particles results from their correlation with 
the reaction plane. Measuring the left-hand side of Eq. (|B|) in various phase space windows, one can then reconstruct 
v n from this equation, up to a global sign. For instance, the E877 Collaboration uses the correlations between three 
rapidity windows to extract flow from their data [fl8f . 

However, nonflow correlations do exist. The two-particle distribution can generally be written as 

/(Pl,P2) = /(Pl)/(P2) + /c(Pl, P2), (6) 

where / c (pi,P2) denotes the correlated part of the distribution. There are various sources of such correlations, 
among which global momentum conservation, resonance decays (in which the decay products are correlated), final 
state Coulomb, strong or quantum interactions ||||]. 

In the coordinate system we have chosen, where the reaction plane is fixed, / c (pi,P2) is typically of order 1/iV 
relative to the uncorrelated part, where N is the total number of particles emitted in the collision. This order of 
magnitude can easily be understood in the case of correlations between decay products, such as p — > mr. A significant 
fraction of the pions produced in a heavy ion collision originate from this decay, and the conservation of energy and 
momentum in the decay gives rise to a large correlation between the reaction products. Since a large number of p 
mesons are produced in a high energy nucleus-nucleus collision, the probability that two arbitrary pions originate from 
the same p is of order 1/N. This 1/N scaling also holds for the correlation due to global momentum conservation [0,0. 

Inserting Eq. (^J) in expression (||), one finds, instead of Eq. (0), 

e in ^-«\ = » n (X>i)v„(X>a) + (e m ^~^) . (7) 
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The left-hand side represents the measured two-particle azimuthal correlation. The first term in the right-hand side 
is the contribution of flow to this correlation, while the second term (e tn (^~My denotes the contribution of the 
correlated part f c . The latter term corresponds to azimuthal correlations which do not arise from flow: we call them 
"direct" correlations, in opposition to the indirect correlations arising from the correlation with the reaction plane, 
that is, from flow. 

Since the correlated two-particle distribution / c (pi,P2) is of order 1/N, so is the second term in the right-hand 
side of Eq. (Q) , which therefore reads 

e «(0i-«\ = v n (V 1 )v n (V 2 )+0(-^). (8) 



t>!xt>2 \N 

However, one must be careful with this order of magnitude. Strictly speaking, it holds only when momenta are 
averaged over a large region of phase space. In the case of the short range correlations due to final state interactions 
(Coulomb, strong, quantum) the correlations vanish as soon as the phase spaces T>\ and T>2 of the two particles are 
widely separated. This is the method used in jy] to get rid of such correlations. If, on the other hand, 2?i and T>2 
coincide, the short range correlations are larger than expected from Eq. (^): in this equation, the total number of 
emitted particles N should be replaced by the number of particles M used in the flow analysis, which is smaller in 
practice. Furthermore, in the case of correlations due to the quantum (HBT) effect, the nonflow correlation scales like 
1 /N only if the source radius R scales like A 1 / 3 |,|. From now on, we shall omit the subscript V for sake of brevity. 
Note, however, that all the averages we shall consider are over a region of phase space which is not necessarily the 
whole space, but may be restricted to the (pT,y) acceptance of a detector. This will be especially important in Sec. 
when we discuss acceptance corrections. 

Equation (||) shows that nonflow correlations can be neglected if v n > N^ 1 / 2 . At SPS energies, the flow is weak and 
this condition is not fulfilled. Indeed, we have shown js]^ that the values of flow measured by the NA49 Collaboration 
at CERN are considerably modified once nonflow correlations are taken into account. 







C. Multiparticle correlations and the cumulant expansion 



The failure of the standard analysis is due to the impossibility to separate the correlated part from the uncorrelated 
part in Eq. (||) at the level of two-particle correlations. The main idea of this paper is to perform this separation 
using multiparticle correlations. The decomposition of the particle distribution into correlated and uncorrelated parts 
in Eq. (|^) can be generalized to an arbitrary number of particles. For instance, the three-particle distribution can be 
decomposed as 

dA 

/(pi,p 2 ,p 3 ) = /c(Pi)/c(P2)/c(P3) 



dpidp 2 dp 3 



+ /c(Pl,P2)/c(P3) + /c(Pl,P3)/c(P2) + /c(P2,P3)/c(Pl) 
+ /c(Pl,P2,P3), (9) 

where / c (Pi) = /(Pi)- The last term / c (Pi, P2, P 3 ) corresponds to the genuine three-particle correlation, which is of 
order 1/N 2 . 

To understand this order of magnitude, let us take a simple example: the to meson decays mostly into three pions. 
First of all, this decay generates direct two-particle correlations: the relative momentum between any two of the 
outgoing pions is restricted by energy and momentum conservation. The corresponding correlation is of order 1/N as 
discussed previously in the case of p — > tttt decays. It corresponds to the second, third and fourth term in the right- 
hand side of Eq. (g). As stated above, the last term in this equation stands for the direct three-particle correlation. 
The corresponding correlation between the decay products of a given u> is of order unity, while the probability that 
three arbitrary pions come from the same ui scales with N like 1/N 2 . Thus the correlation between three random 
pions is of order 1/N 2 . 

More generally, the decomposition of the fc-particle distribution yields a correlated part / c (pi, • • • , Pfe) of order 
1/N k ~ 1 . Generalizing the above discussion of ui — > iririr decay, the decay of a cluster of fc particles will generate 
correlations / c (pi, ■ ■ ■ ,Pfc') with fc' < k. For instance, momentum conservation, which is an effect involving all N 
particles emitted in a collision, produces direct fc-particle correlations for arbitrary fc. 

Such a decomposition is similar to the cluster expansion which is well known in the theory of imperfect gases [fl9| . 
In the language of probability theory, this is known as the cumulant expansion p0[ | . Equations (^) and (|^) can be 
represented diagrammatically by Figs, [l] and|^. In these figures, correlated distributions f c are represented by enclosed 
sets of points, i.e. they correspond to connected diagrams. 
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FIG. 1. Decomposition of the two-particle distribution into uncorrelated and correlated components. The second term in 
the right-hand side is smaller than the first by a factor of order 1 /N. 




FIG. 2. Decomposition of the three-particle distribution. The last term in the right-hand side is of order 1/N 2 relative to 
the first, while the three remaining terms are of relative order 1/N. 



More generally, in order to decompose the fc-point function /(pi, • • • ,Pfc), one first takes all possible partitions of 
the set of points {pi, • • • , Pfc}. To each subset of points {pi 1; • • • , Pi m }, one associates the corresponding correlated 
function / c (piX) " " j Pim)- The contribution of a given partition is the product of the contributions of each subset. 
Finally, /(pi, • ■ • , Pfe) is the sum of the contributions of all partitions. 

The equations expressing the fc-point functions / in terms of the correlated functions f c can be inverted order by 
order, so as to isolate the term of smallest magnitude: 

/c(Pl) = /(Pi) 
/c(Pl,P2) = /(Pl,P2) - /(Pl)/(P2) 
/c(Pl,P2,P3) = /(Pl,P2,P3) - /(Pl,P2)/(P3) ~ /(Pl,P3)/(P2) ~ /(P2, P3)/(Pl) + 2/(pi)/(p 2 )/(p 3 )- (10) 

The cumulant expansion has been used previously in high energy physics to characterize multiparticle correlations: 
it has been applied to correlations in rapidity |2jj and to Bose- Einstein quantum correlations [^,^3|. In these studies, 
the interest was mainly in short range correlations. The use of higher order cumulants was therefore limited by 
statistics: the probability that three or more particles are very close in phase space is small. In this paper, we are 
interested in collective flow, which by definition produces a long range correlation, so that the limitation due to 



statistics is not so drastic. It will indeed be shown in Sec. Ill D that cumulants up to order 6 can be measured, 
depending on the event multiplicity and available statistics. 

We shall deal with multiparticle azimuthal correlations, which generalize the two-particle azimuthal correlations in 
Eq. and can be decomposed in the same way. Referring to the diagrammatic representation in Figs. ^ and ^ we 
shall name the contribution of / c (pi, • • • ,Pfc) to an azimuthal correlation, i.e. the genuine fc-particle correlation, the 
"connected part" of the correlation or, equivalcntly, the "direct" fc-particle correlation. 



D. Measuring flow with multiparticle azimuthal correlations 



Our method, which we now explain, allows the detection of small deviations from an isotropic distribution. If the 
source is isotropic, there is no flow, and the orientation of the reaction plane does not influence the particle distribution. 
We can therefore consider that the reaction plane has a fixed direction in the laboratory coordinate system, so that 
the cumulant expansion can be performed in that frame: in other terms, we replace cj> by the measured azimuthal 
angle <fi. One then measures the fc th cumulant of the multiparticle azimuthal correlation, which is of order N 1 ~ k if 
the distribution is isotropic. Flow will appear as a deviation from this expected behaviour. 

Let us be more explicit. We are dealing with azimuthal correlations. When the source is isotropic, that is, if the 
fc-particle distribution remains unchanged when all azimuthal angles are shifted by the same quantity a, the flow 
coefficients (^) obviously vanish. Therefore, the two-particle azimuthal correlation (J^) reduces to its connected part, 
of order 1/N. As a further consequence of isotropy, averages like /e m ^ 1+ ^ 2 ~^ 3 ^ vanish: only 2fc-particle azimuthal 
correlations involving fc powers of e m ^ and fc powers of e~ m< ^ are nonvanishing. For instance, the fou r-pa rticle 



correlation ^ e m ('t>i+4>2 <fo 4>a)\ j s a priori non vanishing. Introducing the cumulant expansion defined in Sec. II C , this 
correlation can be decomposed into 

^ e in(<t> x +<j>2-<j>3-<l>4)^ — ^ e «"(0i-03)^ ^gin(02-<M^ _|_ ^ e ™(0i-<M^ ^g"l(02-</>3)^ _|_ ^ e i«(0i+02-03-</>4) ^ 



Note that most terms in the cumulant expansion disappear as a consequence of isotropy. The first two terms in the 
right-hand side of Eq. ([ll]) are products of direct two-particle correlations, and are therefore of order 1/N 2 , while the 
last term, which corresponds to the direct four-particle correlation, is much smaller, of order 1/N 3 . However, in the 
case of short range correlations, it may rather be of ord er 1/M 3 , where M is the number of particles used in the flow 
analysis, for the same reasons as discussed in Sec. 



II B 



We name this latter term the "cumulant" to order 4 and denote it by ^ e ™(<i>i+<h-<h-<i>*)y Using Eq. (|l]), it can 
be expressed as a function of the measured two- and four-particle azimuthal correlations: 

e in(0i+02-03-04)\\ = /gi«(0l+02-03-04)\ _ /gi™Oi-<fe)\ /gi«(<fe-04)\ _ / e m(0i-<£ 4 ) \ / e in(02-</>3) \ ^2) 



The reason why we introduce a new notation here is that the cumulant to order 4 will always be defined by (|1^ 
in this paper, even when the source is not isotropic. Now, if the source is not isotropic, the decomposition of the 
four-particle azimuthal correlation involves many terms which have been omitted in Eq. (O) (see Appendix A 1 ) , so 
that the ciimiila.nt // P in (<t>i+4>2-4>3-<t>4)\ nr , l n np-er mrresnonds tn the connected na.rt / P M<h~+<t>2-<p 3 ~<t>4)\ _ 



that the cumulant ^"(tf'i+^-fo-W j n0 i onger corresponds to the connected part (eM^rffo-fo-W^ ^ 

In the isotropic case, the cumulant ^e ln ^ 1+ ^ 2 ~^ 3 ~^y involves only direct four-particle correlations: the two- 
particle correlations have been eliminated in the subtraction. In order to illustrate this statement, let us consider two 
decays p — > irir, and "turn off" all other sources of azimuthal correlations. We label 1 and 2 the pions emitted by the 
first resonance, 3 and 4 the pions emitted by the second. There are correlations between tti and 7T2, or between 7T3 
and 7T4, so that the measured four-particle correlation, i.e. the left-hand side of Eq. (JTl]) , does not vanish. However, 
there is no direct four-particle correlation between the four outgoing pions, so that the cumulant ([l2|) vanishes. More 
generally, if particles are produced in clusters of k particles, there are measured azimuthal correlations to all orders, 
but the cumulants to order k' > k vanish. 

Let us now consider small deviations from isotropy, i.e. weak flow. The two-particle azimuthal correlation receives 
a contribution v 2 according to Eq. (Q). For similar reasons , th e four-particle correlation gets a contribution v*. The 



cumulant defined by Eq. ( |12|) thus becomes (see Appendix A 1 ) 



e*<*+*»-*-**>)) = -vt + (± + V ^j (13) 

where the coefficient —1 in front of v n is found by replacing each factor e m ^ or e _m * in the left-hand side with its 
average value v n . The flow v n can thus be obtained, up to a sign, from the measured two- and four-particle azimuthal 
correlations, with a better accuracy than when using only two-particle correlations, as we shall see shortly. 

It should be noticed that the cumulant involves a contribution from the higher order harmonic 2n, of magnitude 
v 2 n /N 2 . This contribution does not interfere with the measurement of v n provided the following condition is satisfied: 

\v 2 n\^Nv 2 n . (14) 

Since v n is measurable only if v n 3> 1/N, as we shall see later in this section, the interference with the harmonic 
2n occurs only if \vm\ 3> \v n \. In practice, the only situation where this might be a problem is when measuring the 
directed flow V\ at ultrarelativistic energies, where elliptic flow i>2 is expected to be larger than v\. On the other 
hand, this interference will not endanger the measurement of v%, since va should be much smaller. 

In the following, we shall always assume that condition is fulfilled. Then, using Eq. (|l3|), it becomes possible 
to measure the flow v n as soon as it is much larger than i\T^/ 4 . The sensitivity is better than with the traditional 
methods using two-particle correlations which, as we have seen, require v n 3> N~ x / 2 . 

Similarly, using 2fc-particle azimuthal correlations and taking the cumulant, i.e. isolating the connected part (which 
amounts to getting rid of nonflow correlations of orders less than 2k), one obtains a quantity which is of magnitude 
N 1 ~ 2k for an isotropic source. Flow gives a contribution of magnitude u„ . The contribution of higher order harmonics 
Vkn can be neglected as soon as 



\Vkn 



« N k ~ l v k n . (15) 



If \v n \ ^> 1/N, this is not a problem, unless \vk n \ 3> \v n \- This is unlikely to occur, since one expects v n to decrease 
rapidly with n. Neglecting higher order harmonics, there remains the contributions of flow, of magnitude v 2 ^ , and of 
direct 2A:-particle correlations, of magnitude N 1 ^ 211 . Therefore, 2fc-particle azimuthal correlations allow measurements 
of v n if it is larger than N~ 1+1 / 2k . Since k is arbitrarily large, one can ideally measure v n down to values of order 
1/N, instead of 1/yN with the standard methods. A necessary condition for the flow analysis is therefore 

Vn » ^ (16) 



G 



which will be assumed throughout this paper. As we shall see in Sec. [II D, the sensitivity is in fact limited experi- 
mentally by statistical errors due to the finite number of events. 

In practice, the cumulants of multiparticle azimuthal correlations will be extracted from moments of the distribution 
of the Qn-vector introduced in next section. 



III. INTEGRATED FLOW 



In this section, we show how it is possible to measure the value of v n integrated over a phase space region. This 
measurement will serve as a reference when we perform more detailed measurements of azimuthal anisotropies, in Sec. 
IV. We first define in Sec. Ill A a simple version of the Q n -vector, or event flow vecto r, whi ch is used in the standard 
flow analysis to estimate the orientation of the reaction plane. We then show, in Sec. Ill B . that the integrated value 
of the flow can be obtained from the moments of the Q n distribution: eliminating nonflow correlations up to order 
2k by means of a cumulant expansion, we obtain an accuracy on the integrated v n of magnitude jV _1+1 / 2fc , better 
than the accuracy of standar d met hods if k > 1. Instead of using a single event vector Q n , one can do a similar 
analysis using subevents (Sec. Ill Cp . Since the order 2k of the calculation is arbitrary we obtain with either method 
an infinite set of equations to determine v n . The order 2k which should be chosen when analyzing experimental data 
depends on the number of events available (Sec. IIID[ ). More g eneral forms of the Q„-vector, which allow an optimal 
flow analysis, are discussed in Sec. HI E . Finally, in Sec. 1IIF , we recover, as a limiting case, the results obtained in 
the limit of large multiplicity where the distribution of Q n is Gaussian [ ^3|j24| ] . In the whole section, we assume the 
analysis is performed using a perfectly isotropic detector; corrections to this ideal case will be dealt with in Sec. [v|. 



A. The Q-vector 



1. Definition 



Consider a collision in which M particles are detected with azimuthal angles <p\ , • • • , <pM ■ In order to detect possible 
anisotropies of the (f> distribution, it is natural to construct an observable which involves all the </>j, i.e. a global 
quantity. For the study of the n th harmonic, one uses the n th transverse event flow vector ||, which we write as a 
complex number 



Qn 



1 



<M 



M 
3 = 1 



(17) 



where cf>j denotes the azimuthal angle of the j particle with respect to the reaction plane. 

For simplicity, we have associated a unit weight with each partic le in Eq. (|l7|). The generalization of our results to 
arbitrary weights is straightforward and will be given in Sec. HIE. The Q n -vector generalizes to arbitrary harmonics 
the transverse momentum transfer introduced by Danielewicz and Odyniec [g], which corresponds to n = 1 and the 
transverse sphericity tensor introduced in pi| , [25| , which corresponds to the case n — 2. 

In practice, the number of particles M used for the flow analysis is not equal to the total multiplicity ./V of particles 
produced in the collision, since all particles are not detected. However, M should be taken as large as possible. In 
this paper, we shall assume that M and N are of the same order of magnitude. The factor \j\f~M in fron t of Eq. 
(|l7j), which does not appear in previous definitions of the flow vector [0,f3|, will be explained in Sec. Ill A3. 



2. Flow versus nonflow contributions 

A nonvanishing value for the average value of the flow vector, (Q n ), signals collective flow. Indeed, using Eqs. 
and (|l7|), it is related to the Fourier coefficient v n — (e m< ^) by 

(Qn) = VMv n . (18) 

Note that (Q n ) is real, as is v ni due to the symmetry with respect to the reaction plane. 

As stated before, the purpose of the flow analysis is to measure v n , i.e. (Q n )- This is not a trivial task because the 
azimuth of the reaction plane is unknown, so that the phase of Q n is unknown. The only measurable quantity is \Q n \, 
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the length of Q,- 
angles: 



Its square Q n Q* n , where Q* n denotes the complex conjugate, only depends on relative azimuthal 



1 



M 



j-<t>k) 



(19) 



In Sec. Ill B| , we shall see that the flow can be deduced from the moments of the distribution of |Q n | 2 , i.e. from the 
average values (\Q n \ 2k ), where k is a positive integer. To illustrate how flow enters these expressions, we discuss here 
the second order moment (|Qn| 2 ). Averaging Eq. ( |l9| ) over many events and using Eq. (Q), one obtains 



(\Qn\ 2 ) = £ 



M + M(M - 1) (vl + (e in( ^~ 4 ' k) ^ 



(20) 



The first term corresponds to the diagonal terms j — fc, i.e. to "autocorrelations" . If there are no azimuthal correlations 
(neither flow nor nonflow), only this term remains and the average value of \Q n \ 2 is exactly 1. The second term 
corresponds to j ^ k, i.e. to the two-particle azimuthal correlations discussed in Sec. [IB. Since ^e m ^ 1_< ^ 2 ') is 
at most of order 1/M, direct correlations give a contribution which is a priori of the same order of magnitude as 
autocorrelations, although it may be smaller in practice. Equation (pfl) can thus be written 



(\Qn\ 2 ) =M 

= (Qn) 



M \M 
1 + 0(1). 



(21) 



As expected from the discussion of Sec. II B, since (|Q n | 2 ) involves two-particle correlations, flow measurements based 
° n (|Qn| 2 ) are reliable only if \v n \ ^> \j\fM. Smaller values of flow can be obtained using higher moments of the 
distribution of |Q„| 2 , as explained in Sec. [II B. 



If flow is strong enough, the event flow vector can be used to estimate the orientation of the reaction plane. Indeed, 
if |^n| ~> l/y/M, Eqs. Jl§| ) and ( pil| ) show that Q n ~ (Q n ) — \fMv n - Then the phase of Q n is approximately if 
v n > and 7r if v n < 0. Experimentally, one defines Q n as in (|i"7|), with the azimuthal angles <f>j measured with respect 
to a fixed direction in the laboratory (rather than the reaction plane, which is unknown). Then the azimuthal angle of 
the reaction plane <pn can be estimated from the phase of Q n , which we write n(j)Q: 4>R — (resp. 4>r ~ 4>Q + ^ /n) 
modulo 2ir/n if v n > (resp. v n < 0). 



3. Varying the centrality 

Let us now explain the factor 1/y/M in the definition (0). This factor was introduced independently by 
A. Poskanzcr and S. Voloshin, and in |26[ . It is important when using events with different multiplicities M in 
the flow analysis, i.e. events with different ccntralitics. This is the case in practice: one takes all events in a given 
centrality interval in order to increase the available statistics. 

If there is no flow, Eq. ( ^l| ) shows that (|Q„| 2 ) is independent of M since nonflow correlations scale like 1/M. This 
can be understood simply: the sum in Eq. (fi~7]) is a random walk of M unit steps, therefore it has a length of order 
\/M, which cancels out with the factor 1/i/M in front. Flow, on the other hand, depends strongly on centrality (it 
vanishes for central and very peripheral collisions): according to Eq. (|2l|), it gives a positive contribution to (|Q n | 2 ) 
which strongly depends on M. This allows to disentangle flow and nonflow effects. 

Note that flow can be detected by studying the variation of (\Q n \ 2 ) with centrality. This is the method used 
in p7| : one expects (|Q„| 2 ) to be minimum for the most peripheral collisions where the density of particles is too 
small for collective behaviour to set in, and for central collisions where v n also vanishes from azimuthal symmetry. 
However, such a method does not allow an accurate measurement of flow: it is impossible to select true (i.e. with 
6 = 0) central collisions experimentally, and there may still be some flow up to large impact parameters, as suggested 
by hydrodynamic calculations in the case of elliptic flow pij , and by recent measurements P8| . 

The method presented in this paper is more powerful in the sense that it allows flow measurements for a given 
centrality. The error on the centrality selection (due to the fact that one always selects events within a finite range of 
impact parameters) is compensated by the factor \j\[M in the definition of Q n . 
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B. Cumulants of the distribution of \Q n 2 

For sake of brevity, we now drop the subscript n and set n — 1 until the end of this paper, unless otherwide stated. 
All our results can be easily generalized to the study of higher order v n 's by multiplying all azimuthal a ngles by n. 

The moments of the \Q\ 2 distribution involve the multiparticle azimuthal correlations discussed in Sec. \l D . While 
(|Q| 2 ) involves two-particle azimuthal correlations, as seen in Eq. (|l9|), the higher moments (|Q| 2fc ) involve 2fc-particle 
correlations. For instance, we have 



<l^| 4 > = X72 E (e^+^-^A. (22) 



M 2 

These highe r order azimuthal correlations can be used to eliminate nonflow correlations order by order, as explained 
in Sec. II D. This will be achieved by taking the cumulants of the distribution of \Q\ 2 , which we shall soon define. 



1. Isotropic source 



Following the procedure outlined in Sec. |ID|, we first consider an isotropic source (no flow). Using Eq. (|2I 
(|Q| 2 ) is then of order unity, and so are the higher order mome nts ( |Q| 2fc ). However, by analogy with the cumulant 



decomposition of multiparticle distributions introduced in Sec. II C , we can construct specific combinations of the 
moments, namely the cumulants of the Q distribution, which are much smaller than unity: the cumulant ((|Q| 2fe ^ to 
order k, built with the (\Q\ 2 -') where j < k, is of magnitude 1/M . 

As an illustration, let us construct the fourth order cumulant (|Q| 4 )- If the multiplicity M is large, most of the 
terms in Eq. (^2|) are nondiagonal, i.e. they correspond to values of j, k, I and m all different. Then, using the 
cumulant of the four-particle azimuthal correlation defined by Eq. ( |l2] ) and summing over (j, k, I, m), it is natural to 
define (\Q\ 4 \ as 



«M 4 )) = <M 4 )-2<|Q| 2 > : 



(23) 



The order of magnitude of (|Q| 4 )) is easy to derive: each term of type (|l2] ) is of order 1/M 3 as discussed in Sec. |T 
there are M 4 such terms in the sum (p2|); taking into account the factor 1/M 2 in front of the sum, (|Q| 4 )) is finally 
of order 1/M. As intended, two-particle nonflow correlations, which are of order unity, have been eliminated in the 
subtraction (p3|). 

A more careful analysis must take into account diag onal terms for which two (or more) indices among (j,k,l,m) 
are equal. This analysis is presented in Appendix A 2, where we show that diagonal terms are also of order 1/M: 
they give a contribution of the same order of magnitude as direct four-particle correlations. In the following, we 
shall assume that this property, namely that the contribution of diagonal terms is at most of the magnitude of the 
contribution of nondiagonal terms, also holds for higher order moments. 

Among these diagonal terms are the autocorrelations already encountered in the expansion of \Q\ 2 [see the discussion 
below Eq. (|2C|)1, which we define as the term s which remain in the absence of flow and direct correlations. A 

-1/M. 



straightforward calculation (see Appendix A 2) shows that their contribution to the cumulant (|Q| 4 )) is —1/M. As 
in the case of the second order moment (|Q| 2 ) discussed previously, autocorrelations are a priori of the same order 
of magnitude as other nonflow correlations. As we shall see later in this section, they can easily be calculated and 
removed order by order. 




FIG. 3. Decomposition of (|<2| 4 ) = (QQQ*Q*). In the right-hand side, the first term is of order unity while the second term 
is of order 1/M. 



Arbitrary moments (|Q| 2fe ) can be decomposed into cumulants, which can then be isolated in a similar way. This 
deco mpo sition can be represented in terms of diagrams, like the decomposition of the multiparticle distribution in 
Sec. IIP . This is explained in detail in Appendix g. For example, the decomposition of (|Q| 4 ) is displayed in Fig. ||. 







In these diagrams, each dot on the left (resp. on the right) of the dashed line represents a power of Q (resp. Q*), and 
correlated parts, which correspond to direct correlations, are circled: the equation displayed in Fig. [3| stands for 

<|Q| 4 )=2<|Q| 2 )) 2 + <|Q| 4 )). (24) 



Since (|Q| 2 ) = (|<3| 2 ), one recovers Eq. (|23|). More generally, to decompose (IQI 2 *), one draws k dots on each side of 
the dashed line. The diagrams combine all possible subsets of the dots on the left with subsets of the dots on the right 
containing the same number of elements. The latter condition is due to the fact that the average value of (Q l Q* m ) 
vanishes when I ^ m, as a consequence of isotropy. 

In order to invert these relations, and to express the cumulants as a function of the measured moments, the simplest 



way consists in using the formalism of generating functions, recalled in Appendix B2. There, it is shown that the 
cumulant (|Q| 2fe ) is obtained from the expansion in power series of x of the following generating equation, and then 



the identification of the coefficients of 



2k. 



00 2k I 00 2k 



= In (70(2x1(51)) , (25) 



where Iq is the modified Bcssel function of order 0. Expanding this equation to order x 4 , one recovers Eq. (p3|); to 
order x 6 , one obtains the sixth order cumulant 

«M 6 )) = (IQI 6 ) - 9<|Q| 4 ) (|Q| 2 ) + 12(|Q| 2 ) 3 , (26) 
which is of the order of 1/M for an isotropic source. 



2. Contribution of flow 



Let us now consider small deviations from isotropy. As explained in Sec. [ID, these deviations will contribute to 



the cumulants (|Q| 2fc ) defined above. 

The contribution of flow to the fourth order cumulant (|Q| 4 )) is calculated in detail in Appendix [A]. It is shown in 
particular that the diagonal terms in Eq. ( p2] ) are at most of the same magnitude as nondiagonal terms, as in the case 
of an isotropic source. As in Sec. II D, higher order harmonics can be neglected as soon as condition ( pT[ ) is fulfilled. 
One then obtains 



(27) 



where the term — 1/M is the contribution of autocorrelations, i.e. the case j = k = I = m. From Eqs. ( |l8| ) and 
one can measure values of the integrated flow v down to A/~ 3 / 4 , instead of M -1 ' 2 with tradition al m ethods. 

Increased sensitivity can be attained using higher order cumulants. As shown in Appendix B3, the cumulants 
defined by Eq. 



5|) are related to the flow by the following generating equation: 



E 

k=0 



■2k 



(kl) 2 



21, 



}=\nl {2x (Q))+Mhil 



2x 



(28) 



Expanding this equation up to order x 2k , and isolating the coefficient of x 2k /(kl) 2 , one obtains a relation with on 
the left-hand side the cumulant (|Q| 2fe ), while the first term on the right-hand side is the contribution of flow, and 
the second term corresponds to autocorrelations. This identity holds within an error of order M 1 ~ k due to direct 
2/c-particle correlatio ns. U sing Eq. (|l8|), it therefore allows measurements of v within 0(M~ 1+1 / 2fe ), as expected from 
the discussion of Sec. [ID. Expanding Eq. (E8h to order x 4 . 



one recovers 



(p7|). To order x 6 , one obtains 



«IQ| 6 )=4 ( Q) 6 + ^ 



O 



1 

M 2 



(29) 



which extends the limit of detectability down to v ~ M~ 5 / 6 . 

Since Eqs. ( p5| ) and ( p8| ) can be expanded to any order, one obtains an infinite set of equations to determine the 
same quantity (Q). The best choice for the order k will be discussed below in Sec. 
point, we shall discuss an alternative method to measure (Q), the so-called "subevent" method. 



[II D. Before we come to this 
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C. Subevents 



The standard flow analysis, instead of studying the autocorrelation of the event flow vector as in Sec. [II B, deals 
with "subevents" : the set of detected particles is divided randomly into two subsets I and II of equal multiplicities, and 
the two corresponding (subevent) flow vectors Q\ and are constructed. Then one studies the azimuthal correlation 
between Q\ and Qn )2]|| . This is usually done under the assumption that the only azimuthal correlation between the 
subevents is due to flow. Then, from the flow of two equivalent subevents, one can deduce the flow of the whole event 
by a simple multiplication by a factor of y/2, as will soon be explained. 

A nice feature of that method is that, since the subevents have no particle in common, autocorrelations are 
automatically removed: only correlations due to flow and direct correlations remain. Therefore, one may prefer to 
work with subevents when direct correlations are small (although they are, generally, of the same order of magnitude 
as autocorrelations). 

In this section, we shall improve the standard subevent method, in the spirit of Sec 



IIIB 



we eliminate nonflow 

azimuthal correlations order by order by means of a cumulant expansion of the distribution of Q\ and Qn, thereby 
increasing the sensitivity of the method. 



1. Definitions 



Consider two separate subevents of multiplicity M\ and respectively (in practice, one chooses M\ = Me). We 
can construct the subevent flow vectors Q\ n and Qn n as following: 



Mi 



<2lln = 



Mi 
1 



/Mr 



3 = 1 

fc=i 



i s-\ i in^Tt 



(30) 



(31) 



As in Sec. IIIB , we set n = 1 and drop the subscript n in Q\ n and Qn„; generalization to higher n is straightforward. 
By analogy with Eq. (hq), we write 



(Qi) = VMiv h (Q n ) = VMrVr, 



(32) 



where Vj and vr denote the values of v\ associated with each subevent. Hereafter, we shall assume that the two 
subevents are equivalent, i.e. v\ — i>n = v, as is the case if they are chosen randomly. 

Therefore, if M\ — Mn = M/2, the value of (Q) for the whole event is related to the value for one subevent, (Qi), 

by 



(Q) = VMw = V2 (Qi) . 
The purpose here is to measure (Qi), which is equivalent to measuring (Q). 



(33) 



2. Limitations of the standard method 

In order to extract the azimuthal correlation between the subevents, the simplest possibility is to form the product 

Q,Ql = - 7 =L= V = IQiQnl^* 1 -**). (34) 

VMiMn j k 

Using Eq. (Q), the average over many events gives 

{QiQD - VMiMn (v 2 + (eMi-t*)}^ . (35) 



This equation is analogous to Eq. (20), with the important difference that autocorrelations [the first term in the right- 
hand side of Eq. (|2C|)] no longer appear. As a consequence, (QiQJ) vanishes if there are no azimuthal correlations 
between particles. 
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However, two-particle nonflow correlations do remain. Since they are of order 1/N, Eq. ( |35| ) can be written 



{QiQD = VMiMk 

One recognizes in the right-hand side of this equation the flow and nonflow contributions to two-particle azimuthal 
correlations, as in Eq. (||). The only difference lies in the global multiplicative factor \fM\M-&. In particular, summing 
over many particles does not decrease the relative weight of nonflow correlations, as might be believed: they add up 
in the same way as the correlations due to flow. 



O 



N 



(36) 



3. Beyond the standard method 



Now, following the procedure outlined in Sec. Ill B , it is possible to eliminate nonflow correlations between Q\ 
and order by order. This is done by means of a cumulant expansion, which is a trivial generalization of the one 
presented previously. The equation to an arbitrary order 2k is obtained by replacing, in Eqs. (|25|) and (Eq), \Q\ 2 with 



QiQ\, and (Q) with (Qi) (Qe). For example, Eqs. (E3J) and (E7f) become 



(QjQi 2 ) - 2 {QiQlf = MiM n 



0{ IF 



(37) 



which allows measurements of the flow when v is much larger than 1/iV 3 / 4 : the sensitivity is better than with Eq. ( |3(: 
where v is to be compared with 1/iV 1 / 2 . The term — 1/M in Eq. (Bjj), which reflects autocorrelations, is automatically 
removed in Eq. (|37|). 

To make a long story short, the same techniques apply to subevents as to the whole event. The only interest of 
subevents is that they remove autocorrelations. However, they do not remove direct nonflow correlations, which may 
be of the same order of magnitude. Further more, autocorrelations can be also be subtracted systematically when 
working with the whole event, as shown in Sec. IIIB. Another drawback of the subevent method is that each subevent 
contains at most half of the total multiplicity, resulting in increased errors. As a conclusion, the subevent method 
seems to be obsolete when working with cumulants. 



D. Statistical errors 



The cumulant expansion allows in principle the measurement of v down to values of 1 /N by going to large orders 



k, as explained in Sec. II D. In practice, however, since the number of events A cvts used in the analysis is finite, the 
sensitivity is limited by statistical errors. In this section, we determine, as a function of M and N evts , which order of 
the cumulant expansion should be chosen so as to obtain the most accurate value of the integrated flow. 

First, there is a "systematic" error, which is the error due to nonflow correlations. Expanding Eq. ( pq ) to order 2k, 
we obtain an equation relating the measured cumulant (|Q| 2fe )) and the integrated flow (Q), which is of the type 

(\Q\ 2k ) = a k {Q) 2k + 0(M 1 - k ), (38) 

where ak is a numerical coefficient of order unity, and the last term is the systematic error. The resulting error on 
(Q) is therefore 

5(Q) sys t - (Q) 1 - 2k M 1 -" (39) 

The systematic error thus decreases with increasing k, since (Q)yfM = Mv n 3> 1, as assumed in Eq. (|l6|). 

Let us now discuss the statistical error. When averaging a quantity over a large number of events N cvts , the 
statistical error is generally of relative order 1 / y/ N cvts . Since the moments of the distribution of \Q\ 2 are of order 
unity, the absolute statistical error on the moments is of order l/y/N cv t s . The same error applies to the cumulants, 
which are constructed from the moments. If there is no flow, a more accurate calculation shows that the statistical 
error on the cumulant is 

V-^evts 
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If the flow is weak, that is if (Q) <C 1, this formula still holds approximately. Using Eq. (paj), one thus derives the 
statistical error on the integrated flow 

S(Q)st,t-(Q) 1 - 2k N- v 1 t / s 2 . (41) 

Since we have assumed that (Q) <§; 1, the statistical error increases with k, unlike the systematic error. 

It is very likely that this property still holds in the more general case when (Q) is not much smaller than unity. 
However, we have not been able to derive a general formula for the statistical error for arbitrary (Q) and k. We only 
have formulas for the lowest order cumulants. Using the cumulant to order 2 (k = 1), 



1 1 + 2(Q) 2 

o (Q) stat = tttTa \ — 77 , (42) 



2(Q> V 



and with the fourth order cumulant (k = 2), 



*<QW = ^- • (43) 



V CVtS ' 



One sees on these two formulas that for very strong flow ((Q) 3> 1), the statistical error 8(Q) sta t is of order 1 / \J N e . 
independent of (Q). This remains true for higher order cumulants. Note, moreover, that both formulas give Eq. ( f4l| ) 
in the limit of small (Q). 

Since the systematic error decreases with k and the statistical error increases with k, the best accuracy is achieved 



for the value of k such that both are of the same order of magnitude. Using Eqs. ( |39| ) and (41), one thus obtains the 
optimal value of the order 2k: 

Since, in practice, M is at least of the order of a hundred at ultrarelativistic energies, the fourth order cumulant 
2k = 4 [i.e. Eq. (p7|)] gives the best accuracy if the number of events lies in the range 10 3 < iV cv t s < 10 6 . Higher 
order cumulants may be useful if a large statistics is available and/or if the multiplicity M is low, as for instance in 
a peripheral collision. 

The flow is detectable only if (Q) is larger than both statistical and systematic errors. Taking for instance N cvts = 
10 5 and M = 300, statistical and systematic errors are of the same order. One then obtains, using Eq. (f43|), that 
flow can be seen if (Q) > 0.3. Using Eq. (|lj), v can be measured down to 1.6% using the fourth order cumulant. If 
V = 3%, a typical value at the CERN SPS, then (Q) ~ 0.5. Using Eq. @, the typical error is then 5{Q) ~ 0.02, i.e. 
6v = 0.1%. 



E. Weighted Q-vectors 

The vector Q n has been defined in Eq. ( |l7| ) with unit weights. A more general definition is 



M 



1 Vw.-e*"*', (45) 



where the weight Wj is an arbitrary function of pt, y, the particle type, and the order of the harmonic under study. 
As a consequence, we shall restore the index n in this subsection. 



1. Flow analysis with arbitrary weights 



The method exposed in Sec. IIIB also applies with this more general definition. There are only two slight differences. 
The first is that the average value of Q n , which we have denoted by (Q n ), is no longer related to the average value v n of 
the flow by Eq. ([l8]). This modification is not important for what follows: we shall see in Sec. IV B that measurements 
of differential flow depend on the value of (Q n ) rather than v n . The second difference is that autocorrelations cannot 
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be removed so simply: the procedure given in Appendix B 4 is no longer valid, so that the subevent method, which 

avoids autocorrelations, may regain some interest. 

Apart from this difference, the procedure is the same as in Sec. IIIB. In particular, the cumulants of the event flow 
vector distribution are expressed in the same way in terms of the moments. The generating equation (|2£ 



still holds, 

with the caveat that the last term, corresponding to autocorrelations, is no longer exact. 

However, autocorrelations are unchanged at the lowest order: a calculation analogous to the one leading to Eq. (|2(i| ) 
shows that (|<3n| 2 ) = 1 if there are no azimuthal correlations between particles, up to terms of order 1/M. Changes 
occur only at higher orders. 



2. Optimal weights 

What is the best choice for the weight w{pT,y,n)l In practice, it should be chosen so as to maximize the effect 
of flow: one should try to obtain a value of (Q n ) as large as possible, since this value will determine the accuracy 



in the measurement of azimuthal distributions, as we shall see in Sec. IV. From the definition ([15|), averaging over 



azimuthal angles and denoting by (v n )j the value of v n for the corresponding particle, one obtains 



, . _ £ 3 -=iK) 3 -w 3 - 



M 

^(Vn)l (46) 



where we have used a simple triangular inequality, and the fact that the flow coefficients (v n )j are real. The identity 
holds when Wj — X(v n )j, where A is arbitrary. In other terms, the optimal weight for a particle with given rapidity 
and transverse momentum is the associated flow coefficient (v n )j itself. 

Of course, since the goal is precisely to measure v n , the above discussion does not answer the question of the choice 
of the optimal weight. However, general properties of the u n 's can be used to guess a reasonable choice of w. Since 
v n is an odd (resp. even) function of the center of mass rapidity for odd n (resp. even n), so should be w. Regarding 
the pt dependence, one may note that at low pt, v n generally behaves as v n cx pVj^ Q. Therefore, it seems natural 
to choose w cx p^, when measuring the n th harmonic. For n = 1, Q n then becomes the sum of transverse momenta, 
weighted by an odd function of rapidity, which was the definition chosen in Q. For n — 2, Q n is then equivalent to 
the transverse momentum sphericity tensor used in p4|j. 



F. Gaussian limit 

In this section, we compare our method to methods previously used in jl3|,^4]], which rely on the large multiplicity, 
Gaussian limit. It is well known that, according to the central limit theorem, the distribution of the fluctuations of 
Q around its average value (Q) is Gaussian in the limit of large M. Up to corrections of order 1/M, the normalized 
probability of Q = Q x + iQ y thus reads 



d 2 Q 2ira x fJ v \ 2a x 



with a 2 x = (Ql) - (Q) and a 2 y = (Q 2 y ). 

We shall first show that this limit is equivalent to the cumulant expansion to order 4 presented in Sec. IIIB. Then 
we shall discuss the relationship with an alternative method to measure flow, which has been used in the literature, 
and consists in fitting the distribution of \Q\. 



1. Higher harmonics 



In the case of the Gaussian distribution (J47|), one easily calculates the cumulants used in Sec. Ill B 



(|Q| 4 ) - 2 (|Q| 2 ) 2 = - <Q) 4 + 2(a 2 x - a 2 ) <Q) 2 + (a 2 - a 2 ) 2 



(48) 
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In order to compare these equations with Eqs. (J2l]) and (A7), we need to evaluate the sum a 2 = cr| + a 2 and the 
difference a 2 — o 2 . 



From Eqs. Jig) and ( |20| ), one obtains 

^ 2 = (IOI 2 )-(Q> 2 

= 1 -i>? + (M- 1) (e^-**) 



(49) 



where the last term is of order unity since <^e* 
vector (ft5[). One thus recovers Eq. 
Let us now calculate the difference: 



) is of order 1/N < 1/M. This still holds for the generalized 



1 M 

°"x - = T7 H K cos (<& + <W) _ (COS <j>j) (C0S( 



M 

v 2 - v\ + 0(v 2 ), 



(50) 



where the first two terms in the last equation come from the diagonal terms j = k, while th e re maining term is the 
contribution of nondiagonal terms. Reporting this expression into Eq. (fig), we recover Eq. (A7): higher harmonics 
reflect a deviation from isotropy in the fluctuations of Q. 



2. Isotropic fluctuations 



Neglecting higher harmonics, we may write a x = a y . Then the distribution ( ft7| ) becomes 



dp 

d^g 



■ exp 



\Q - (Q) 



(51) 



With this distribution, we expect to recover the results of Sec. IIIB, where higher harmonics were also neglected. 
Indeed, one finds after some algebra, for arbitrary, real x, 



ln(I Q (2x\Q\)) = <r 2 x 2 +\nl (2x (Q», 



(52) 



to be compared with Eqs. J25| ) an d (ffiij ). According to Eq. J49|), the extra term a 2 x 2 is of order unity, in agreement 
with the statement following Eq. (|28|) that the correction at order x 2k is 0(M 1 ~ k ). 

Corrections to the central limit theorem are of order 1/M. Thus, expanding Eq. ([32]) in powers of x, one obtains 
identities which are valid up to that order. To order x 4 , we recover the result obtained in the previous section, see 
Eq. (p7|), with the same accuracy. To order x 2k with k > 2, the results obtained in Sec. [II B 
we have seen that the correction is of magnitude M 1 ~ k <C 1/M. 



are more accurate since 



3. Distribution of \Q\ 

A method for extracting the flow from the data, which was proposed in JD||24|], consists in plotting the measured 
distribution of \Q\. This method led to the first observation of collective flow in ultrarelativistic nucleus- nucleus 
collisions [Ml]. It is a simplified version of the method based on the sphericity tensor p9[, which led to the first 



observation of collective flow at Bevalac |30|. Note that these methods are more reliable than what we call the 
"standard method" in this paper, in the sense that one need not neglect nonflow correlations. 
The distribution of |Q| is obtained by integration of Eq. (|l]) over the phase of Q: 




1 dp 

One must then fit both parameters a and (Q) to the data. 

If there is no flow, that is (Q) = 0, the \Q\ distribution given by Eq. (|5|) is purely Gaussian: 



(53) 



1 dp 2 / 

■exp . (54) 



T 2 
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The \Q\ distribution deviates from the Gaussian shape if the flow is strong enough compared to the fluctuation scale, 
that is for values of (Q) > a. In particular, the maximum of the distribution is shifted to \Q\ ^ if (Q) > a. Since 
a is of order 1, using Eq. (jl^), this condition is equivalent to v > \ j\f~M. Note, however, that one need not assume 
v 1/vM, as with the methods based on two-particle azimuthal correlations. 

If (Q) <C a, i.e. v -C 1/yM, the shape of the distribution is very close to a pure Gaussian distribution. In fact, 
the deviations from the Gaussian shape are of order (Q) /ct 4 [pi. This can be seen by expanding Eq. ([53]) to order 
(Q) 2 , which is equivalent to replacing a 2 with a 2 + (Q) 2 in Eq. (|54|). Alternatively, one can eliminate a and obtain 
(Q) directly using the following identity, which can be easily derived from Eq. j5l]): 

(\Q\ i )-2(\Q\ 2 ) 2 = -{Q)\ (55) 

again showing that the deviation is of fourth order in the flow. Knowing that the deviation to the central limit is of 
order 1/M, one finds that Eq. (|55|) is equivalent to Eqs. (^3|) and ( p7j), i .e., to the cumulant expansion to order 4. 

The importance of the factor l/\/M in the definition of Q, Eq. (|17|), also appears clearly when fitting Eq. ( |5^ ) 
to experimental data. Because of this factor, a does not depend on the multiplicity M in the limit of large M, as 
discussed above. This is especially important when the fit is done using events with different multiplicities M. If 
there is no flow, the distribution of |Q| is Gaussian with width a. If a depended on M, the distribution would rather 
be a superposition of Gaussian distributions with different widths. In this case, the left-hand side of Eq. ( |55| ) would 
be positive, hiding a possible weak flow. This phenomenon probably explains why the first analysis of the E877 
Collaboration juj gives zero values of the flow in some centrality bins. 

When fitting Eq. ( |53| ) to the data, it is important to fit independently (Q), which reflects the flow, and c, which 
also involves two-particle correlations, according to Eq. (^9|). Assuming that a is the same for all Fourier harmonics, 
as was done by E877 |Q, amounts to neglecting two-particle correlations. 

Finally, note that the Gaussian limit can also be applied to the subevent method, yielding interesting results: 
in particular, the distribution of the relative angle between Q\ and is not the same for direct correlations and 
correlations due to flow |2(J . 



IV. DIFFERENTIAL FLOW 



In this section, we explain how it is possible to perform detailed measurements of azimuthal distributions: typically, 
one wishes to measure v n for a given type of particle as a function of the rapidity y and the transverse momentum 
Pt- In the following, we shall call this particle a "proton," but it can be anything else. We denote by ip its azimuthal 
angle, and by v' m the corresponding differential flow coefficients v' m = (e mi,/i ). Unlike the standard method, as stated 
before, we do not make the assumption that all azimuthal correlations are due to flow. As in the case of the integrated 
flow studied in Sec. Ill, we get rid of nonflow correl ations order by orde r, by means of a cumulant expansion. 

The principle of the method is explained in Sec. [V A . In Sec. [VB , we show that v' m can be obtained from the 
azimuthal correlation between ijj and the flow vector Q. As in the case of integrated flow, the order to which n onflow 
correlations must be eliminated depends in practice on the number of events available: this is explained in Sec. IV C , 
where we also estimate the resulting accuracy on v' m . Our method is compared to traditional methods in Sec. IV D. 



A. Principle and orders of magnitude 

The differential flow coefficients v' m can be obtained only through azimuthal correlations with other particles, 
typically particles used to estimate the orientation of the reaction plane, which we call "pions" in this section, 
although they can be anything else. 

For instance, correlating the proton with one pion, v[ can be derived from the measurement of the two-particle 
azimuthal correlation 

(eW-^)=v[v 1 + o(±y (56) 

where v\ refers to the pion, and is determined independently. We have used an analogy with Eq. (||). The term 0(1/N) 
comes from two-particle nonflow correlations between the proton and the pion. The error made in the determination 
of v[ is thus of order l/(Nvi). Of course, one should correlate the proton to particles with a strong flow, so that vi 
be as large as possible. 



16 



More accurate measurements can be obtained using higher order correlations and a cumulant expansion. For 
instance, at fourth order, one can eliminate the two-particle nonflow correlation by correlating the proton with three 
pions and taking the cumulant, by analogy with Eqs. ( |l2] ) and (13): 



e i{ip+4>l-4>2-4>3) \\ = / e i(V>+0l-02-03) \ _ / e i(il>-<t>l) \ / e «(0i-03)\ _ / e i(ip-4>3)\ / e «(0l-02) 



More generally correlating the proton with 2k + 1 pions, the connected part of the correlation is of order 1/N 2k+1 
[since it corresponds to direct (2k + 2)-particle correlations], while the contribution of flow is v[v 2k+1 . Comparing 
both terms, the accuracy on v[ is thus of order 1/ (Nvi) 2k+1 . Using Eq. (|l^), this shows that the accuracy increases 
with increasing fc, i.e. when using multiparticle correlations. 

Higher harmonics, such as v' 2 , can be obtained by at least two methods. The first consists in multiplying all the 
angles by 2 in the equations above, and replacing v[ and v\ with v' 2 and v 2 , respectively. A second method is to 
mix two different harmonics, measuring source is isotropic, this quantity is of order 1/N 2 since 

it involves a direct three-particle correlation. If there is flow, neglecting other sources of correlation for simplicity, 
^(2V-0i-02)^ f ac torizes into (e 2 ^) (e - ^ 1 ) (e - ^ 2 ) = v' 2 v\. Putting everything together, we obtain 

eW-^W)=^ + o(l). (58) 

One sees that nonflow correlations come into play only at order 1/N 2 , rather than 1/iV when comparing the same 
harmonics as in (|5^). Nonetheless, they do not disa ppear . These correlations can also be eliminated order by order 



using the cumulant expansion, as we shall see in Sec. IV B . Generally, if one correlates the proton with 2k + m pions, 
one obtains an accuracy on v' m of order l/(Nvi) 2k+m . 



B. Differential flow from correlations with Q n 



In order to correlate a proton with pions, it is convenient to use the event flow vector Q„, Eq. ([l7|). From now 
on in this section, we choose n = 1, and drop the subscript n, i.e. we write Q and v instead of Q\ and V\. On the 
other hand, we keep the subscript m for the proton v' m because several harmonics may be measured. Generalization 
to arbitrary n is straightforward: one simply multiplies all azimuthal angles (of both protons and pions) by n. 

In the standard flow analysis, one usually excludes "autocorrelations" by excluding the "proton" under study from 
the definition of the event flow vector that is, the azimuthal angle ip is not one of the <pj in Eq. ([l7|). Within our 
method, one can still do so, but it is not even necessary. First, aut ocorr elations will be removed order by order as well 
as direct correlations, as in the case of the integrated flow in Sec. IIIB . Furthermore, autocorrelations, if any, can be 
subtracted e xact ly if the event flow vector Q is defined with unit weight, as in Eq. (|l7j). This subtraction is performed 
in Appendix C4. For simplicity we neglect the corresponding term in this section, unless otherwise specified. 

Let us start with the measurement of the first harmonic v[. The two-particle azimuthal correlation between the 
proton and a pion, Eq. (p6|), can be expressed introducing the vector Q defined by Eq. (|l7]). Summing Eq. ( |56| ) over 
all the pions involved in Q, one obtains the correlation between (Q) and the proton: 



(QV*) = (Q) 



O 



1 

Nv 



(59) 



The value of (Q) must be obtained independently using the methods discussed in Sec. III. 

More accurate measurements, involving correlations of the proton with several pions, are performed using higher 
order moments, as in Sec. IIIB. These higher order moments are obtained by weighting the previous expression with 



powers of |Q| 2 , i.e. by measuring (\ l Q\ 2k Q*e l ^'Y These moments are then decomposed into cumulants. For instance, 
Eq. ( p7| ) becomes 



<IQI 2 Q ; 
-{Qf 



o 



2(( 

\(N V y 



)(M 2 ) 



(60) 



through which we define the cumulant ( ! \Q\ 2 Q*e 1 ^^) . 

As in the case of integrated flow, the decomposition of higher order moments {\Q\ 2k Q* e 1 "^) in cumulants can be 
represented in terms of diagrams. For instance, the decomposition of (\Q\ 2 Q*e 1 ^) is displayed in Fig. ||. 
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FIG. 4. Decomposition of (\Q\ 2 Q*e 1 ^} = (Q* 2 Qe 1 ^ . The cross correspond to the factor e 1 ^ of the proton, while dots on 
left (resp. on the right) of the dashed line stand for Q (resp. Q"). 



The diagrams in this figure stand for 



(|Q| 2 Q*e^) = 2 <Q*e^» «|Q| 2 )) + {\Q\ 2 Q*e^} 
= 2<QV^)<|Q| 2 ) + «|Q| 2 QV^> 



(61) 



One thus recovers the expression of the cumulant, Eq. (|6(j). More generally, in order to decompose the moment 
(\Q\ 2k Q*e 1 ^) — (Q k Q* k+1 e 1 ^), one draws a cross on the left representing the proton, k dots on the left and k + 1 dots 
on the right representing the pions. The graphs combine all possible subsets of the points on the left with subsets of 
the points on the right containing the same number of elements. 

Let us now discuss the measurements of higher harmonics of the proton azimuthal distribution v' m . In the case 
rn = 2, Eq. (|58|) gives, summing over the pions involved in Q, 



(Q* 2 e 2 ^) = (QY 



O 



( 1 



(62) 



To obtain a better accuracy, one must decompose higher order moments ^|Q| 2fe Q* 2 e 21 ^) in cumulants. In terms of 
the diagrammatic representation, the proton is now associated with two crosses, as seen in Fig. ^]for k = 1. 



x 
x 



= 3 



FIG. 5. Expansion of (jQ| 2 Q* V"^) = (Q* 3 Qe 2 ^). The linked crosses stand for the proton, while dots on left (resp. on the 
right) correspond to Q (resp. Q*). 



As before, the graphs combine all possible subsets of the points on the left with subsets of the points on the right 
containing the same number of elements, with the subsidiary condition that the two crosses belong to the same subset. 
In the left-hand side of Fig. ||, the dot on the left of the dashed line can be associated with any of the three dots on 
the right. The equation represented by the figure can be written 



(|Q| 2 Q*V^) =3«Q* 2 e 2 ^))«|Q| 2 )) + ( 

= 3(Q* 2 e 2 ^)(\Q\ 2 ) + (\Q\ 2 Q* 2 e 2 ^} 



(63) 



where the last term involves a direct five-particle correlation, and is therefore of order M 2 x 0{\/N A ). When there 
is flow, one obtains 



2 Q* 2 e 2 »M 



(IQI 

(Q) 4 



*2 2i^p 



)-3(Q* 2 e 2 ^){\Q\ 2 ) 



~2Vn + O 



(Nv) 4 



(64) 



Cumulants of arbitrary order, for arbitrary harmonics v' m , can be obtained by expanding in powers of x the following 
generating equation, derived in Appendix C2: 



E 



„2fc+m 



fc! (k + m)\ 



(l m (2x\Q\)(Q*/\Q\) m e^) 
(Io(2x\Q\)) 



(65) 



where I m is the modified Bessel function of order m. For m = 1, one recovers Eq. ( |6l| ) by expanding this equation to 
order x 3 . For m = 2, one recovers Eq. (|63|) by expanding this equation to order x 4 . 
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The cumulants defined by Eq. (|6^) are related to the differential flow by 



00 „2fe+m 

EX 
, _ o k\ (k + m)\ 



M 



To (2x/ 



(66) 



The second term corresponds to autocorrelations, and must be included only if the proton is involved in the flow 
vector Q. In the case m = 1, one recovers the lowest order formulas (59) and ( |60| ) by expanding this equation to order 
x and x 3 , respectively. For m = 2, one recovers Eqs. 



At order x 2k+m , Eq. (pq ) gives an accuracy on v' ri 



EVA. 



]62| ) and (|64|) by expanding it to order x 2 and x 4 , respectively, 
of order l/(Nv) 2k+m , as expected from the discussion of Sec. 



C. Statistical errors 

Equation (^) generates an infinite set of equations to measure the differential flow v' m , since it can be expanded 
to any arbitrary order x 2k+m . As in the case of integrated flow, the best choice of k is the one that yields the best 
accuracy on v' m . It results from a compromise between systematic errors stemming from nonflow correlations, which 
decrease when using higher order cumulants, and statistical errors, which increase with the order k. 

The equation obtained when expanding Eq. ( |66| ) to order x 2k is of the type 

(\Q\ 2k Q* m e im ^) = b k (Q) 2k+m v' m + O (M-*-( m / 2 )) , (67) 

where bk is a numerical coefficient of order unity. Neglecting for the moment the error on the integrated flow (Q), 
this equation gives a systematic error on v' m 

(Sv' m ) syst ~ (Q)- 2k - m M- fe -(™/ 2 >. (68) 

This systematic error decreases when increasing the order k. Note that for m = 1, the systematic error should be the 
same on the differential flow as on its integrated value at a given order. Thus we expect Eqs. (p8|) and (^) to give 
the same result, using (|l8|). In doing the comparison, one must pay attention to the fact that the cumulant used for 



differential flow \\Q\ Q*e i v\ involves 2k + 2 particles while the cumulant for integrated flow « Q 2fc ]> involves only 
2k particles. Thus, comparing the two "at a given order" means that we must replace 2k in Eq. (38) by 2k + 2 in Eq. 
©■ 

The statistical error on the cumulant ( |67| ) is of order 1 / y/N{, vts , where N{, vts is the number of events containing a 
proton. This leads to an error on v' m : 

(Sv'J stat ~ (Qy 2k - m (N' vts )- 1/2 (69) 

If m = 1, we again recover the result obtained for the integrated flow, provided we replace 2k by 2k + 2 in Eq. 
and N^ vta by MN cvts (which is the total number of particles involved in the measurement of integrated flow) in Eq. 

If (Q) < 1, the statistical error ( |69|) increases with increasing k, and the optimal value of k is that for which 
statistical and systematic errors are equivalent, i.e. 

2fc~ -m+ lnN Z ts - (70) 
In M y ' 



The error on (Q), estimated in Sec. [II D, should also be taken into account. However, the measurement of differential 
flow is done in a limited region of phase space, by definition, so that the corresponding statistics is smaller than for 
the integrated flow where many more events can be used. It is then safe to assume that the statistical error on (Q) 
gives a negligible contribution to the error on v' m . 

If m = 1, the previous equation shows that k = 1 is more accurate than k — (the latter value corresponds to 
the standard method, neglecting correlations) only if lnN^ vts / InM > 2, i.e. if the statistics is large enough, typically 
iVg Vts > 10 4 for an event multiplicity M ~ 100. For higher harmonics m > 1, the contribution of nonflow correlations 
are smaller as explained above: thus the lowest order method k = is to be chosen unless a very large number of 
events is available, typically iV cvts > 10 6 for the second harmonic m = 2 if M ~ 100. 
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D. Relation with previous methods 



Previously used methods fiplj also study the correlation between t he eve nt flow vector (Q7h with the momentum 



of the proton. The traditional justification is that, as explained in Sec. Ill A, the phase n(j)Q of the event flow vector 
( |l7| ) gives an estimate of the orientation of the reaction plane modulo 2ir/n. Studying the correlation between ip and 
4>q, one can reconstruct harmonics v' n , v' 2n , v' 3n , etc. 

The standard analysis relies on a purely angular correlation. One measures the average (cosm(^ — 4>q)). Neglecting 
nonflow correlations, this quantity is the product of v' m and a resolution factor which is given by an independent 
measurement |32| . Our method relics on similar averages, weighted by powers of \Q\: 

(|Q| 2fe+m cosm(V> - 4>q)) = (Q* m+k Q k e tmip ) . (71) 

In the traditional method, autocorrelations are usually removed explicitly by specifying that the proton under study 
is not used in constructing Q n in Eq. ([l?]) g. However, nonflow direct correlations, which are of the same order of 
magnitude as autocorrelations, do remain, and limit the sensitivity of the analysis. With our method, autocorrelations 
can be removed in the same way as in the standard analysis. But we also remove direct correlations, thereby increasing 
the sensitivity of the measurements. 

V. ACCEPTANCE CORRECTIONS 

For simplicity, the discussion has been limited so far to an ideal detector, i.e. a detector with an acceptance which 
is azimuthally isotropic in (f>. An actual detector is never perfect, either because its components are of uneven quality, 
or simply because it does not cover the whole 4> range. In this section, we discuss a simple extension of the method 
which allows us to work with any detector. More precisely, it allows the detection of deviations from an isotropic 
source, i.e. flow, with any detector, and the correction is implemented in the same way for all detectors. However, the 
accuracy on the measurement of v n can be poor if the detector covers only a limited range in tj>. 



The only modification lies in the definition of the cumulants, for which the expressions given in Sees. Ill B and 



IV B are no longer valid. These modified cumulants are defined in Sec. VA for integrated flow and in Sec. VB for 
differential flow. As we shall see, the analytical expression of higher order cumulants become very lengthy, so that it 
is more convenient to work directly at the level of gen erating functions. As an illustration of our method, results of a 



simple Monte-Carlo simulation are given in Sec. V C 



A. Integrated flow 

The key idea is that anisotropies in the detector acceptance can be handled much in the same way as anisotropies 
of the emitting source. The only difference is that the relevant coordinate system is the laboratory system in the first 
case, and the system associated with the reaction plane in the second case. 

Let us be more specific: until now, we have been working in the coordinate system associated with the reaction 
plane, i.e. with the emitting source. In this system, we used a cluster exp ansio n to define direct fc-particle correlations, 



of order A^ 1 k relative to the uncorrelated fc-particle distribution (Sec. |HC|). This cluster expansion allowed us to 
construct the "connected moments" of the distribution of Q, which were noted as (Q k Q* l ) c (Appendix B 1), of order 
M 1 relative to the corresponding moment {Q k Q* ls j. This decomposition was performed for an arbitrary source, 
but with an ideal detector. 

Exchanging the roles played by the source and the detector, the same reasoning applies if we work with an isotropic 
source and an imperfect detector, provided we use the coordinate system associated with the detector. We thus define 



the connected moments exactly in the same way, replacing by the measured (see Sec. II A). Similarly, the flow 
vector Q will be denoted Q when azimuthal angles are measured in the laboratory system, i.e. when cj>j is replaced 
by <f>j in the definition (|l7j). 

I f the acceptance is not perfect, averages such as (e m ^) or (e™^ 1 ^ 2 ^ 3 )) no longer vanish. Thus, nondiagonal 
moments (Q k Q* 1 ) with k ^ I are also nonvanishing: there is no more cancellation due to isotropy, and all terms must 
be kept in the cumulant expansion. At order 2, for instance, the cumulants are defined as 

m ^ (q 2 ) - (q)\ 

))^(\Q\ 2 )-(Q)(Q*). (72) 
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These cumulants are of the same magnitude as when the acceptance is perfect, i.e. of order unity, while the moments 
(Q 2 ) and (|Q| 2 ) scale like the multiplicity M if the detector is very bad. Note that at this order (k + I = 2), taking 
the cumulant is equivalent to shifting the distribution of Q by its average value (Q), as proposed in ||. 

Higher order cumulants can be obtained in similar way as for an ideal detect or. The only difference is that the 
simplifications due to isotropy no longer exist. Thus one cannot use expression (pa) for the generating function of 



the moments; one must use instead the more general expression (B2). The cumulant (Q Q* ) is therefore defined by 



E 



kj 



Z Z 



(Q k Q* 1 ) = InGo(z) = ln(e» ) . (73) 



Expanding the right-hand side to order z* k z l , one obtains the cumulant as a function of the measured 

moments (Q k Q* 1 ) with fc' < fc and V < I. While the moment (Q k Q* 1 ) is of magnitude M^ k+l ^ 2 for a bad detector, 
the corresponding cumulant (Q k Q* 1 ) is of order M(k+i)/2 N i-k~i _ M i-(fc+0/2_ 

If the acceptance is not too bad, we assume that relation ( |28| ) between the cumulants and the integrated flow is 
approximately preserved. The integrated flow can then obtained from the cumulants to order 2, 4, 6 using Eqs. ( pT 
( P7| ) and ((2^), which we write again in the form: 




(Q) =(IQI 2 I))- i + ±\ " ■ ( 74a ) 



(0) ^- (I , --L + (-L) ± y 1 + W^ + W , 

(«>* = i wo -i^GK^vfcr (74c) 

where, in the right-hand side of each equation, the last three terms stand for autocorrelations, systematic errors due 



to direct 2fc-particle correlations, and statistical errors due to the finite number of events (see Sec. [II D ), respectively. 
Note that (Q) denotes the average value of Q in the coordinate system associated with the reaction plane, i.e. what 
we call the "integrated flow". It must not be mistaken for (Q) [see for instance Eq. (|7^)1, which denotes the average 
value in the laboratory coordinate system, and vanishes if the acceptance is perfect. Note also that only the "diagonal 
cumulants" (|0| 2fc )) (i- e - with fc = I) are related to the flow. These diagonal cumulants could equivalently be written 
(IQ| 2fc ) since Q and Q differ only by a phase. Other cumulants, {Q k Q* 1 } with k I, are not influenced by the 
flow and vanish except for statistical and systematic errors. They can therefore be used to estimate the magnitude of 
errors. 

The modified definition of higher order cumulants involves a large number of terms when the detector acceptance 
is nonisotropic. For instance, the fourth-order cumulant is obtained by expanding Eq. ([73]) to order z 2 z* 2 : 

«IQ| 4 ) = (IQI 4 ) - 2 (Q) (QQ* 2 ) - 2 (Q*) (Q*Q 2 ) - 2 (|Q| 2 ) 2 - (Q 2 ) (Q* 2 ) 

+ 8 (Q) (Q*) (\Q\ 2 ) + 2 (Q) 2 (Q* 2 ) + 2 (Q*) 2 (Q 2 ) - 6 (Q) 2 (Q*f . (75) 



This equation replaces Eq. (|23J) for an imperfect detector. It shows that implementing acceptance corrections order 
by order can be very tedious since it involves a large number of terms. 

It is simpler to work directly with generating functions. Although this might seem to be more complicated, it is not 
unnatural since the generating functions constructed from experimental data have the same geometrical properties 
as the data, in particular regarding the detector acceptance. For instance, when the detector is isotropic, so is the 



generating function Eq. (B5) 



One can compute numerically the generating function of the cumulants Go(x,y) at various points in the complex 
plane, then extract numerically the coefficients at a given order by means of an interpolating polynomial. Let us be 
more specific: separating the real and imaginary part of the flow vector, we write it as 



M 



Qx = , = J2 W 3 COS ( n< ^ )> 

'£.,•=! 3=1 



M 



® y = / m E w 3 sin^j ) • ( 76 ) 

T i=1 W 2 j = l 
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The generating function of the cumulants, defined by Eq. (|7^), is a real- valued function: 

]nG (x,y) = l n ^e 2xQ " +2yQ y^ , (77) 

where we have set z = x + iy. 

According to Eq. (73), the cumulant to order 2k, (|Q| 2fc } is the coefficient of (zz*) k = (x 2 + y 2 ) k in the power 
series expansion of this generating function, up to a factor l/(fc!) 2 : 



]nQo(x,y) = ^ 



(78) 



where we have kept only the relevant terms in the expansion. The cumulant can be obtained from the tabulated 
values of hiQo(x,y) using the interpolation formulas given in Appendix D 1. 



B. Differential flow 



When measuring differential flow, acceptance corrections can be implemented in the same way as for integrated 
flow. Flow is extracted using the same formulas as when the detector is perfectly isotropic in azimuth (Sec. |rv| ), 
without the simplifications all owe d by isotropy. T here fore, one must take as the generating function of the cumulants 



C m {z) the general expression (|C4|) instead of Eq. (jC6J) . We thus define the cumulants by 

= z* Q+zQ* +irmp 



£ w {w***)) = m*) - ' (79) 

k,l \ / 

where ip denotes the azimuthal angle of the proton, measured in the laboratory coordinate system. This equation 
replaces Eq. ( |65|) for an imperfect detector. Expanding Eq. (|79| ) to order z for m — 1, we obtain for instance 



- (Q*) (e**) (80) 



We assume that the relation (g^) between the cumulants and the differential flow, v' m , is approximately preserved 
if the acceptance is not too bad. For k — and k = 1, flow is then related to the cumulants by Eqs. ( |59] ) and ( |60| ) for 
m = 1 and by Eqs. (js^) and (Q) for m — 2. We rewrite these formulas: 

(Q) v[ = «QV<*» + O (j^j ± -^=, (81a) 

(Q)\[ = - (|Q| 2 Q*e^» + o (^) ± -^=, (81b) 



( Q ) 2 ^ = <Q* 2 e^» + O (^-J ± - 7 ==, (81c) 
(Q> 4 4 ^ -5 «]Q| 2 Q* 2 e 2 ^» + O (^) ± (81d) 



where, in the right-hand side of each equation, the second term represents the systematic error due to direct particle 
correlations, while the last term is the statistical error due to the finite number of events. Note that only the cumulants 
(Q k Q* ! e""*} with I = k + m are related to the flow by Eqs. (|l|). 

We wish to recall here that the differential flow v' 2 might also have b een o btained from the correlation between the 
azimuth of the proton and the event flow vector Q2. As stated in Sec. [V B| , the only modification is a multiplication 



of all ang les by 2, so that this does not change Eqs. ( |79| ) and j82|). Therefore, v' 2 may be deduced from Eqs. (81a 
and (|81b| ) by the simple substitution of v[ and (Q) by v' 2 and (Q2) respectively. 

As in the case of integrated flow, the modified definitions of the cumulants quickly involve a large number of terms 
when going to higher orders. Therefore, it is simpler in practice to extract the cumulants numerically from the 
generating function. For this purpose, one must tabulate numerically the real and imaginary parts of C m (z): 
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e 2xQ !c +2yQ y C0S ( TO ^;) 

$t[C m (x,y)\ = 



/ e 2xQ x +2yQ y ^ 
' £ 2xQ !c +2yQy sm ( m? £) V 

*lC m ( X ,y)}= * {e2xQx+2yQy) ' ■ (82) 
Keeping only the terms with I = k + m which are related to the flow, the generating function (|79|) becomes 



\2k r)*m £> imip\ 



fc=0 



Interpolation methods to calculate the cumulants ((|Q| 2fc Q * m e tTra ^]) as a function of the tabulated values of the 



generating function, are explained in detail in Appendix D 2 



C. Results of a Monte-Carlo simulation 



We have tested our method with a simple Monte-Carlo simulation. Particles have been generated randomly with 
the distribution 

dN 

— cx 1 + 2vi cos + 2v 2 cos(2</>) . (84) 
d(p 

The value of the integrated directed flow, which we tried to reconstruct, was fixed to v% = 0.03, corresponding roughly 
(up to a sign) to the value measured at SPS for pions JT^. We have taken various values of t>2, in order to probe the 



interference between both harmonics, discussed in Sec. II D 



In a first step, we do not simulate nonflow correlations between the particles. In order to take into account the 
effect of detector inefficiencies, we have assumed that all particles are detected, except in a blind azimuthal sector of 
size a. The simulation has been performed with AT evts = 200000 events, and a multiplicity M = 200 for each event. 
For simplicity, we have assumed that exactly 200 particles are detected in each event. Fluctuations in M should not 



influence the results, as explained in Sec. Ill A 3. With these values, the optimal sensitivity for the integrated flow is 



obtained f or k = 2 according to Eq. (kf4J), i.e. by taking the fourth order cumulant. We therefore reconstruct the flow 



using Eq. (|74b[ ) 



With the values we have chosen, (Q) = vivM ~ 0.42 < 1, so that traditional methods might fail, as stated 
before. Within our method, the statistical error on v±, calculated with Eq. (f43|), is of the order of 0.14%. Since direct 
correlations between particles are not simulated, the only systematic error comes from detector inefficiencies and the 
higher harmonic 

Results are shown in the table below. The table gives the reconstructed Vi as a function of the size of the blind 
angle a, and the higher harmonic 





a = 0° 


a = 45° 


a =90° 


a =135° 


a =180° 


V2 = 


3.04 


3.10 


3.11 


2.91 


2.11 


«2 = 3 


2.83 


2.85 


2.98 


2.78 


2.57 


V2 = 6 


2.65 


2.82 


2.78 


3.55 


4.24 


V2 = —3 


3.30 


3.22 


3.23 


2.99 


2.57 



TABLE I. Results of a Monte-Carlo simulation. All values of vi and V2 are given in %. 
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If t>2 = 0, the reconstructed value is compatible with the theoretical value within statistical errors, except for the 
highest value of a, i.e. when the detector covers only half of the range in azimuth. Therefore, errors due to acceptance 
imperfections are under good control. 

The systematic error from hig her harmonics, on the other hand, is far from negligible. The limits of applicability 
of our method, given by Eq. (AS), are here —0.43 < t>2 < 0.07. We have checked these bounds numerically. The value 
V2 = 0.06 is very close to the upper bound. However, the corresponding relative error on V\ is only 12% with an ideal 
detector. 

In a second step, we simulate nonflow correlations: for simplicity, we do this assuming that particles are emitted 
in pairs, both particles in a pair having exactly the same azimuthal angle. This would be the case for the two-body 
deca y of a very fast resonance. Taking the same numerical values as above, the standard method, corresponding to 
Eq. (74a), gives V\ = 7.7%: it f ails, as expected, overestimating the flow by more than a factor of 2. On the other 
hand, the fourth-order formula ( [74 bj ) , which eliminates two-particle nonflow correlations, gives v\ = 3.1%, in much 
better agreement with the theoretical value. 



VI. SUMMARY 



We have proposed in this paper a new method for the flow analysis, which is more sensitive than traditional methods 
to small anisotropics of the azimuthal distributions. In this section, we summarize the procedure which should be 
followed in practice. 



The first step consists in measuring the "integrated flow," as explained in Sec. III. This corresponds to the problem 
of the reaction plane determination in the standard flow analysis. One first constructs, event by event, the flow vector 
Qn defined by Eq. ([76]), where the (f)j are th e azim uthal angles of the particles in the laboratory coordinate system. 
The weight Wj is chosen as explained in Sec. |lll E 2 ; ideally, it should be taken equal to the differential flow v n (pT, y), 
i.e. proportional to p^, and even (resp. odd) in the rapidity y for even (resp. odd) n. Alternatively, one may choose 
the simpler version with unit weights (0). The value of n depends on the system under study: up to energies of 
10 GeV per nucleon, one usually works with n = 1, i.e. with Qi 18 3l| . At SPS, directed flow is so small that a 
better accuracy is obtained by working directly with the second harmonic, i.e. by constructing Q2 fl7j| . Then, only 
even harmonics can be measured. Most of this paper has been written assuming n = 1. In order to generalize the 
results to the case n = 2, one need only multiply all azimuthal angles by 2. 

Measuring the integrated flow amounts to measuring the average value of the flow vector, (Q n ), in the coordinate 
system where the reaction plane is fixed. The average value (Qn) is of order v n y/M (it is even equal to that value if 
one is working with unit weights), where v n is the Fourier harmonic of order n, and M the number of particles used 
in the flow analysis. As explained in Sec. Ill, the integrated flow (Q n ) is obtained from the cumulant (|Qn| 2fe )), which 
removes nonflow correlations up to order 2k, the standard method corresponding to the lowest order, k = 1. The 
value of k is chosen so as to obtain the best sensitivity. It results from a balance between systematic and statistical 
errors, and depends both on the number of events A^ cv t s available for the flow analysis, and on the number of particles 
used to determine the reaction plane in each event, M. The optimal order k is then given by Eq. (|44|). However, 
performing measurements with other values of k does not cost much and provides a useful comparison. 

The cumulant ((|Qn| 2fc ) is a combination of the moments of the distribution of Q n , i-e. it is expressed as a function 
of the measured moments (Q l n Q n m ), with I < k and m < k. In this paper, we have used the formalism of generating 
functions to derive the corresponding formulas at arbitrary order. As explained in Sec. |v|, this is not only an elegant 
formalism: it is also the simplest way to calculate the cumulants numerically from experimental data. For this 
purpose, one tabulates the generating function Qo(x,y), defined by Eq. ([77]), at various points in the (x,y) plane. 
In this equation, the brackets denote an average over the whole sample of events. The cumulant (|<3n| 2,c )) is then 
obtained by extracti ng n umerically the coefficient in front of (x 2 + y 2 ) k in the power series expansion of h\Qv(x,y), 
as explained in Sec. VA. The integrated flow (Q n ) is finally obtained from the cumulant using Eqs. (|74|). 

The value of (Q n ) is the important parameter in the flow analysis, since it determines the accuracy of the recon- 
struction of azimuthal distributions. If (Q n ) > 1, the flow can easily be studied with traditional methods, although 
the present method should give more accurate results. If (Q n ) < 1, on the other hand, standard methods fail, while 
our method still works. 

The second step in the flow analysis is to perform detailed measurements of the flow coefficient v' m for a particle of 
given rapidity and transverse momentum, i.e. differential flow. The coefficient v' m can be obtained from the comparison 
of the azimuth of the particle under study with an event flow vector, which can be either Q m , calculated with the 
same harmonic, or a Q n , calculated with a different harmonic, provided m i s a m ultiple of n. For instance, v' 2 can 
be measured with respect to Q\ or Q2, as explained in ||. We show in Sec. IV E| that it is the value of (Q n ) which 
determines the accuracy on the measurement of v' m . Therefore, n should be chosen so that (Q n ) be as large as possible. 
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For instance, at RHIC where V2 is expected to be much larger than v%, v' 2 should be measured with Q2 rather than 
with Qx, as is already the case at SPS J!?]]. In the text, we have assumed n = 1. If one uses Q2, then m must be 
replaced by 2m everywhere in our equations. 

As the integrated flow, the differential flow v' m is obtained from a cumulant ^\Q\ 2k Q* m e lm ^y which eliminates 
nonflow correlations up to an arbitrary order 2k + m, the standard analysis corresponding to the case k = 0. Here 
again, the best choice of k is the one which leads to the smallest error: its value is given by Eq. (f70|). In order to 
measure the cumulants, one first tabulates the generating function $$2\ ) at various points in the complex plane. The 
cumulant is then obtained by extracting the coefficient proportional to z * k z k + m ' m the power series expansion of the 
generating function, as explained in Sec. VB. Finally, the differential flow v' m is related to the cumulants by Eqs. (|8l|). 

A limitation of our method at a given order is the possible interplay of higher harmonics in the measurement. 
For instance, Eq. ( |lg| ) shows that in the fourth-order cumulant, the second harmonic «2n interferes with v n . More 
precisely, \v2n\ must be small compared with Mv\ (see Eq. (|l4|)). This limitation means that the method should be 
used with much care when extracting the directed flow (n = 1) at RHIC and LHC p4| , since it is expected to be 
much smaller than elliptic flow. On the other hand, in the case n — 2, there should be no problem since V2 is much 
larger than 114. 

While higher harmonics or statistical errors may limit the use of the method, there is no problem with the acceptance 
of detectors. As a matter of fact, the required corrections appear in a natural way in the method, at all orders, from 
a modification of the generating equation which is the same for all detectors. In particular, the sensitivity remains 
unchanged when acceptance corrections are taken into account, so that choosing the order in the expansion of the 
generating equation does not depend on that problem. 

Most of our results have been established in the limit where azimuthal anisotropies are weak. For this reason, our 
method seems to be more adapted to ultrarelativistic energies, i.e. at SPS energies and beyond, where V\ and i>2 
are usually less than 0.1. In particular, it should be very useful in the forthcoming flow analyses at the Brookhaven 
Relativistic Heavy Ion Collider. 
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APPENDIX A: DETAILED STUDY OF THE FOUR-PARTICLE AZIMUTHAL CORRELATION 



In Sec. Al, we calculate the cumulant of the four-particle azimuthal correlation, introduced in Sec. II D. Then, in 



Sec. k 2, we calculate the fourth order cumulant of the Q distribution, introduced in Sec. [II B 



1. Cumulant of the four-particle correlation 



The cumulant of the four-particle azimuthal distribution has been defined by Eq. (|12j) when the source is isotropic. 
We set n = 1 for simplicity: 



e i(01+02-<fe-<M \\ = / e i(<t>l+4>2~<t>3-4>4) 



e i(<Pl-<p3) \ / e i{4>2-4>4.)\ _ / e i{<Pl-<Pi) \ / e i(4>2-4>3) 



(Al) 



Here, we want to evaluate the right-hand side of this equation when the source is no longer isotropic. 

In order to do so, we expand the four-particle distribution into conn ecte d parts, as explained in Sec. [IC . Using 
the diagrammatic representation introduced there, the quantity in Eq. ( |Al[ ) can be decomposed as in Fig. |(| 
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FIG. 6. Expansion into connected parts of the cumulant of the four-particle azimuthal correlation. Dots on the left (resp. 
right) of the dashed line represent e 1 ^ (resp. e - "^). 

The diagrams in Fig. rJ stand for: 



-*(03 + 04) 



(A2) 



Note that the direct two-particle correlations (< 
connected part of the correlation, i.e. (e 



i(01-03) 



} c are automatically removed. In the isotropic case, only the 



i(<t>l+<p2 — <p3— <§>i) 



) , remains in the right-hand side of Eq. (A2) 



Let us now enumerate the orders of magnitude of the different terms in the right-hand side of Eq. (A2). As 
stated above, all terms but the last vanish in the isotropic case: indeed, {V^ 1+ ^ 2 ~^ 3 )) c and ^e ±J ^ 1+< ^ 2 ^) c are not 
invariant under the transformation <pj — > <pj + a, where a is any angle. Therefore, it seems reasonable to consider that 
these terms are proportional to v\ or V2, depending on whether a factor e ±la or e ±2la appears under the previous 
transfor matio n. Furthermore, since we consider here connected /c-particle correlations, they behave like 0(l/iV fc_1 ) 
[see Sec. II C|. More precisely, 



f>i+4 



5 ±iOi+0 2 ) 



o 



(A3) 



Note that the second term in the right-hand side of Eq. ( A2) is smaller than either the first or the third terms. 

Finally, the order of magnitude of the right-hand side of Eq. (Al) is v\ + 0(v 2 /N 2 + 1/N 3 ). We have neglected 
v 2 /N 2 since it is smaller than either vf or 1/A^ 3 . 



2. Calculation of the cumulant 



In this section, we derive the order of magnitude of the fourth order cumulant of the \Q\ 2 distribution, defined by 
Eq. (p3|). From the definition of the event flow vector Eq. (|l7|), one obtains 

(|Q| 4 )) = ^2 [(e l{ ' p ' + ' pk ~" pl ~" pm) ^ - (V ( ^'" 0i) ^(V (0fc " " l) ^ - ^ e i(0J '~ 0m) ^e i(0iI ~ 0i) ^ . (A4) 

j,k,l,m 

In the above sum, one may distinguish nondiagonal terms, when all four indices are different, and diagonal terms, for 
which at least two indices are equal. 

Nondiagonal terms corres pon d precisely to the cumulant of the four-particle correlation. The corresponding con- 
tribution, evaluated in Sec. |A 1|, must be multiplied by the combinatorial factor M(M — 1)(M — 2) {M — 3) ~ M 4 



With the factor 1/M 2 in front of the sum in Eq. (AA), the contribution of nondiagonal terms to (|Q| 4 )) is of order 
M 2 vf + 0(v 2 + l/N). 

We are now going to show that diagonal terms give a contribution at most of the same order as nondiagonal terms. 
Let us enumerate the various diagonal terms: 

i) If j = k = I = 77i, each term in the sum is equal to — 1. This is the contribution which we call "autocorrelations". 
Multiplying by a combinatorial factor M and by the factor 1/M 2 in Eq. (A4), the corresponding contribution is 
exactly -1/M. 



ii) When three indices are identical while the fourth is different, i.e. in 4M(M — 1) cases, the difference in Eq. ( A4) 
reduces to — (V^ 1 "^ 2 '). Using Eq. (g), this contribution is of order — 4vf + 0(1/N). Although this contribution 
is a two-particle correlation, it is suppressed by the combinatorial factor: vf is much s mall er than the term 



M 2 vf which appears in the cumulant of the four-particle azimuthal correlation (see Sec. Al). Therefore, this 
contribution will be negligible. 

iii) Let us consider the cases when the indices are equal two by two. 

• If j = k and I = m but j ^ I, which occurs M(M — 1) times, the difference is given by 

(jWi-to)^ _ 2 ^i-^y =v 2 + ^(^-03) ^ - 2 {y 2 + (e**-^ y . (A5) 



2G 



The order of magnitude is then v\ + 0(1/N). Here, we have neglected terms of order v\/N and 1/N 2 , 
smaller than 1/N; the term vf is smaller by a combinatorial factor 1/M 2 than the similar contribution of 
nondiagonal terms. Note that the higher harmonic v 2 contributes here. We shall see below that these higher 
harmonics can limit the use of our method. 

• The 2M(M — 1) cases {j = m and k — I but j ^ /} or {j — I and k = m but k ^ 1} yield a contribution 

- (e^ 1- ^) . Its order of magnitude is —2v\ + 0(1/N 2 ), negligible compared to nondiagonal terms. 

iv) There are two cases when three indices are different: 

• If j = I or j = to or k = I or k = to, while the two remaining indices are different, the contribution is 

— (V^ 1- ^ 3 ') , to be multiplied by a combinatorial factor 4M(M— 1)(M — 2). Thus, the order of magnitude 
is M[—4vf + 0(1/N 2 )) and this contribution is suppressed by a factor 1/M with respect to the cumulant of 
the four-particle correlation. 

• If the two identical indices are either (j,k) or (l,m), the combinatorial factor is 2M(M — 1)(M — 2), 

which multiplies a term (e^ 2 * 1 -^ 3- ^ 4 )) - 2 (e*^ 1- * 5 )) 2 . Using Eq. (|), the three-particle correlation 
^ e i(20! -03-04)^ can b e expanded as 

/ ^(2^-03-^4)^ _ ^i^x^g-i&j^g-i^ + ^ e i(20i-0 3 )^ ^ e -#4^ _|_ ^ e i(20i-04)^ ( e -«fe^ 

_|_^ e 2i0^ ^ e ~i(03 + 04)\ ) + ^ e i(20 1 -0 3 -04)\ ) 

_ 1w « + (d) + o(k) +0 (^). C (A e) 

The second term in the difference, —2 (V^ 1 "^ 3 )) 2 , gives a contribution of — 2v* + 0(1/7V 2 ). Finally, since 
terms such as vf, vf/N are suppressed because of the combinatorial factor, the contribution in this case is 
2Mv 2 vf + 0{Mvl/N) + 0(M/N 2 ). 

We shall assume that the total multiplicity in the collision TV and the number M of particles used to calculated the 
flow vector are large and of the same order of magnitude. Then, we find that the contribution of the diagonal terms 
is -1/M + vl + 2Mv\v 2 + 0{vl + 1/M). 

All in all, when we add the contributions of diagonal and nondiagonal terms, we obtain the following result: 

«M 4 )) = -jrr m2 < + 2Mv " V2 + v * + °(^ 2 ) + o (Jj) . (A?) 

The first term in the right-hand side coresponds to autocorrelations, the last two terms are due to nonflow correlations, 
and the three remaining terms arise from flow. One would like —M 2 v\ to be the dominant flow term. However, higher 
harmonics, i.e. V2, also contribute. If v 2 is large enough, it may even reverse the sign of the contribution of flow to 
^|Q| 4 )). This does not happen provided v 2 lies in the following interval : 

- Mv 2 (V2 + 1) < v 2 < Mv 2 (V2 - 1). (A8) 



We have checked these bounds with our Monte-Carlo simulation, see Sec. V C 



APPENDIX B: A GENERATING EQUATION FOR THE INTEGRATED FLOW 



In this Appendix, we first construct the cumulants of the distribution of \Q\, which we note 
of t he measured moments (|<5| 2fc ) (Sees 



2k \ 



Bl and B2) 



B3), and show how to remove autocorrelations at all orders (Sec. B4). 



as a function 

Then, we r elate the cumulants to the integrated flow (Q) (Sec. 



1. Cluster decomposition of the moments 



We have shown in Sec. II C how the fc-particle momentum distribution can be decomposed, in a coordinate frame 
where the reaction plane is fixed, into a sum of terms involving lower order distributions (fc' particles with k' < k), 
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plus a "connected" term of relative order l/N k 1 . This decomposition al so ap plies to the moments of the distribution 
of the event flow vector Q defined by Eq. (0). As pointed out in Sec. [11 B . moments of order k involve fc-particle 



azimuthal correlations. This allows us to write a series of equations similar to Eqs. (||) and (^): 

(Q) = (Q) c ( Bla ) 

(Q 2 ) = (Q)l + (Q 2 ) c (Bfb) 

(QQ*) = (Q) c (Q*)c + (QQ*)c (Bic) 

(Q 3 ) = (Q)l + 3 (Q) c (Q 2 ) c + (Q 3 ) c (Bfd) 

(Q 2 Q*) = (Q)l (Q*) c + 2 (Q) c (QQ*) C + (Q 2 ) c (Q*) c + (Q 2 Q*) C , etc. (Bic) 

In these equations, the subscript c denotes "connected" moments. The connected moment of order k is of magnitude 
M 1_fc / 2 : a factor M 1 "* 5 comes from the fact that it involves direct A;-particle correlations (see Sec. |[I C| ) , and a factor 
M k l 2 from the scaling of Q with the number of particles like \/M, see Eq. (|l7|). 

The expansion of a given moment (Q k Q* lX ) in connected parts can be represented graphically by the expansion of 
a [k + ZVpoint diagram into connected diagrams. This is similar to the decomposition of the fc-particle distribution 
in Figs. [I] and g. To be more specific, the decomposition of ( x Q k Q* 1 ) is represented by drawing k dots of one type 
corresponding to powers of Q and I dots of another type corresponding to powers of Q* . One then takes all possible 
partitions of this set of k + I points. To each subset of points one associates the corresponding connected moment. 
The contribution of a given partition is the product of the contributions of each subset. Finally, (Q k Q* 1 ) is the sum 
of the contributions of all partitions. Figure ^ represents, as an example, the decomposition of (Q 2 Q*)- 



• — 






FIG. 7. Decomposition of (<5 2 Q*) in connected parts, see Eq. (Ble). Dots on the left of the dashed line represent factors of 
1 while dots on the right represent factors of Q* . Circled subsets correspond to connected moments. 



The connected moments can be expressed as a function of the moments by inverting Eqs. (Bl) order by order. 
However, this procedure is very tedious. An elegant and compact way to express moments of arbitrary order in terms 
of the connected parts, and to invert these relations, consists in using generating functions. The generating function 
of the moments is a function of the complex variable z which is defined as 



Go(z) 



_ I e z*Q+zC 



k,l 



(B2) 



where k and / go from to +oo, and the brackets denote an average over many events. It is well known in graph theory 
that the generating function of connected diagrams is the logarithm of the generating function of all diagrams p5[ . 
Therefore, the generating function of the connected moments is the logarithm of the generating function of the 
moments 



k.l 



kill 



(Q k Q* l ) c = lnG (z). 



(B3) 



The normalization coefficient l/k\l\ has been chosen such that (Q k Q* l ) c appears with a unit coefficient in the 
expansion of (Q k Q* 1 ), as in Eq. ( |Bl[) . Expanding Eqs. ( [B2| ) and (B3) to order z* 2 z, one finds for instance 



(Q 2 Q*) C = (Q 2 Q*) - (Q 2 ) (Ql - 2 (Q) (QQ*) + 2 (Q) 2 (Q*) 



(B4) 



which can be checked by inverting Eqs. (Bl) order by order. Note that we are working in a coordinate system w here 
the reaction plane correponds to the x-axis, and is unknown. In this coordinate system, the generating function (B2) 
is not a measurable quantity. 
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2. Isotropic source 



We now consider specifically an isotropic source, i.e. without flow. In that case, the moment (Q k Q* 1 ) vanishes if 
k 7^ I. The connected parts (Q k Q* l ) c enjoy the same property. Therefore, in the diagrammatic expansion, one only 
retains terms containing as many powers of Q as of Q* , i.e. as many dots on the left as on the right. The quantity 
represented in Fig. ^ does not satisfy this property, and therefore it vanishes. A decomposition with nonvanishing 
terms is represented in Fig. ^|. 

Keeping only the terms k — I, the generating function (B2) becomes 



0o CO 



oo I >2k 

Ep5<IOI a *> = < J o(2|«0|)>, 



fc=0 



(B5) 



where Iq is the modified Bessel function of order 0. Note that now the generating function Go itself is isotropic, since 
Go {z ) = Go(ze ta ). The consequence is that it can be evaluated in the laboratory coordinate system rather than in 
the coordinate system associated with the reaction plane: it thus becomes a measurable quantity. We define the 
cumulants through 



E Su2 W) = h ^o(z) = In (W\zQ\)) 



k=l 



(fc!) 



(B6) 



They coincide with the connected moments (|Q| 2fc ) c defined in Eq. (B3) if the source is isotropic. Note that for 
an isotropic system, the raw moment (|Q| 2fe ) is of order unity, as noted in Sec. 1IIB. The corresponding cumulant 
(IQI 2 ' 8 } is of order M 1 ~ k . Eq. (B6) corresponds to Eq. (p5|) , where we have set x = \z\. 



3. Flow 



Let us now calculate the cumulants in the case of collisions with flow. Neglecting for simplicity nonflow correlations 
between particles, we can write (Q k Q* 1 ) = (Q) k (Q*) = (Q) k+l . The generating function (B2) thus becomes 



(B7) 



Now, we want to compare with the experimental value of Go(z), which is measured in the laboratory coordinate system 
wher e th e azimuth of the reaction plane <fin ^ 0. The generating function in this coordinate system is deduced from 
Eq. (B7) by the substitution z — > ze l ^ R . Averaging the new expression over all possible (f>R, under the assumption 
that the distribution of cj>R is uniform, one obtains: 



Go (z) 



1 

2^ 



^** + '''- t **M>d4> R = I (2\z\(Q)). 



Gathering the results obtained in Eqs. ( |B6|) and (|B8|), we obtain: 

(\Q\ 2k )=\nGo(z)=lnI (2\z\(Q)). 



h ( fc! ) 2 



(B8) 



(B9) 



Expanding Eq. (BE) to order |z| 2fc , one obtains an equation relating (Q) to the cumulant (|Q| 2fe )). However, when 



writing Eq. (B7), we have neglected direct 2fc-particle correlations and autocorrelations. As explained in Sec. Bl 



both give a contribution of magnitude M 
corrrection of order M 1 ~ k . 



l-k 



to the cumulant (|Q| 2fc ). Thus, Eq. (B9) at order \z\ 2k is valid up to a 



4. Removing autocorrelations 



Equation (B9) can be somewhat refined. In the case of a Q vector defined with unit weights, as in Eq. 
autocorrelations can be calculated and subtracted explicitly, which is the purpose of this section. 
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This calculation has already been done in Sec. Ill for the lowest orders k = 1 and k — 2: we have seen in Eq. (EC 



that diagonal terms give a contribution 1 in the expansion of (|Q| 2 ). In this paper, we refer to these diagonal terms 
as "autocorrelations". Similarly, they give a contribution — 1/M to the fourth order cumulant (|Q| 4 ), see Eq. (|27j ) 
and Appendix |A|. 

To calculate the contribution of autocorrelations to the cumulant at an arbitrary order, we once again make use of 
the generating function Qo{z), Eq- (B2). Neglecting correlations for simplicity, the contributions of the M particles 
to Qo{z) factorize, leading to 



Go(z) 



,(2x cos 4>+2y sin <p) / \/~m\ 



M 



(BIO) 



where we have set z — x + iy, and the brackets here denote an average over <j). Assuming for simplicity that the 
distribution is isotropic, one obtains 



Go(z) = 



-i M 



(Bll) 



This is the expression of the generating function if there are only autocorrelations (no direct correlations, no flow). 
If there is flow, we assume that autocorrelations and flow give additive contributions to the cumulants, which yields 
instead of Eq. @: 



E 

fc=0 



7 \2k f nl 

-' (|Q| 2fe ) =lne (z) = ln/ (2N(Q})+MlnV 



(ft!) 



(B12) 



This formula is equivalent to Eq. (p8|), which we use in Sec. |HIB , It removes exactly all autocorrelations when the 
event flow vector Q is defined with unit weights, as in Eq. (|17|). 



APPENDIX C: A GENERATING EQUATION FOR DIFFERENTIAL FLOW 



In this ap pendix, we follow closely the same procedure as in Appendix |^, applied to differential flow. In Sees. 
CI and C2, we first construct the relevant cumulants (|Q| 2fc Q*'e"™^)), as a function of the measured moments 
(Q k Q* l e lm ^'Y Here, ip denotes the azimuthal angle of the particle under study (which we call a proton), and m the 
order of the harmonic measured for this particle . T hen, we relate the cumulants to the integrated flow v' m (Sec. 
and show how to remove autocorrelations (Sec. C4). 



1. Cluster decomposition 



A quantity such as (Q k Q* l e lm ^ involves correlations between k + I + 1 particles: k + I "pions" (according to 
the terminology introduced in Sec. IV) and a proton. This quantity can therefore be decomposed, in the coordinate 
system where the reaction plane is fixed, into a sum of terms involving lower order correlations, plus a connected term 
of relative order l/N k+l . For instance, we can write 



(Cla) 
(Clb) 

(QQV™") = (Q) c (Q*) c (e im ^) c + (QQ*) C (e»^) c + (Q) c (Q*e im ^) c + (Q*) c (Qe^') c + (QQ*e m ^) c , (Clc) 



(Qe im ^) = (Q) c (e lm ^) c + (Qe lm ^) 



where, in the third equation, the last term is of order 1/iV 2 relative to the first one. Such decompositions can 
be represented diagrammatically, in a way similar to the decomposition of (Q k Q* 1 ) in Appendix We choose to 
represent the proton by m crosses on the left, f or re asons which will become clear below, when we consider the specific 
case of an isotropic source. For instance, Eq. flClq 1 ) can be represented diagrammatically by Fig. @. 



• = 
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FIG. 8. Decomposition of (QQ* e 21 ^^ in connected parts, see Eq. (Clc). As in Fig. [?], the dot on the left (resp. right) of the 
dashed line stands for Q (resp. Q*). The linked crosses represent the proton, the number of crosses being chosen equal to the 
harmonic under study, here m = 2. Circled subsets (connected diagrams) correspond to connected moments. 



In order to express in a compact way the relations between the moments (Q k Q* l e tm ^^ and the corresponding 
connected moments (Cj k Q* l e lm ^) , we introduce the following generating function 



Gm (z) 



}+zQ imip 



v—* Z* k Z l 



^ kin 

k,l 



(C2) 



Expanding Q m (z) to order z* z , one obtains all the moments (Q k Q* 1 e 1 ™"^) . In order to obtain the generating function 
of the connected moments, we note that each diagram in Fig. |8| can be written as the product of a connected diagram 
containing the crosses, i.e. the proton, times an arbitrary diagramfnot necessarily connected) invo lving o nly p ions, 
which correspon ds to the terms (Q k Q* 1 ) considered in Appendix [b|. For instance, using Eqs. ( Bla ) and ( Blc ), one 
can rewrite Eq. (Clc) as 



(QQV"^) = (QQ*) (e m ^) c + (Q) (QV"^) c + (Q*) (Qe m ^) c + (QQ*e m * p ) c 



(C3) 



Therefore, the generating function of the diagrams with pions and protons (e z Q +Z Q \ j s the product of the 
generating function of graphs with only pions, i.e. Go(z) defined in Eq. (B2), by the generating function of connected 
graphs with pions and protons. This latter is therefore: 



C m {z) 



sr^ z* k z l 



^ kill 

k,l 



(Q k Q* l e mip ) 



Qo(z) 



/ e z*Q+zQ*^ 



(C4) 



As in Eq. (B3), the normalization coefficient l/k\l\ has been chosen so that (Q k Q* l e lm ^) appears with a unit 



coefficient in the expansion of (Q k Q' 



2. Isotropic source 



We now consider the particular case of an isotropic source, without flow. Then the moment (Q k Q* l e lm ^) vanishes 
when k + m ^ I, and so do the corresponding connected parts. This is the reason why we chose to represent the 
proton with m crosses: in the isotropic case, only diagrams with the same number of points (crosses and dots) on 
each side of the dashed line do not vanish. 



Expanding the generating function (C2) and keeping only the nonvanishing terms, one finds 

\z\ 2k z m 



Grn(z) = ^ 



fc=0 



k\ (k + m) 



■{\Q\ 2k Q 



) 



I m (2\zQ\) 



\zQ\ 



(C5) 



where I m is the modified Bessel function of order m 

We define the cumulants ( 
source is isotropic. Using Eq 



2kQ*m e i m ip\ g0 tfiak they coincide with the connected moments in Eq. (C4) when the 



this gives 

oo 

C m {z) = ^2 

fe=0 
In 



\z\ 2k z 



(I 



2k Q*m i m ip\ 



kl (k + m)\ 

t (2|*QI)(a) m e' m *) 
W\zQ\)) 

This equation is equivalent to Eq. (p5|), setting x = \z\. 



(C6) 
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3. Flow 



Finally, we turn to the more general ca se o f collisions with flow. Neglecting for simplicity nonflow correlations 
between particles, the generating function (|C2|) becomes 



(C7) 



As explained in Sec. B3, this quantity is measured in the laboratory coordinate system, therefore one must replace z 



by ze**" and average the new expression over all possible That yields 



2^ 
I m {2\z\ (Q)) 



z 

N 

Using Eq. (pcf), the generating function of cumulants ( |C4| ) takes the form 

7 m (2|z|(g)) ( z 



I (2\z\(Q)) \\z 



Gathering Eqs. (C£) and (p9), we obtain 



E 

fe=0 



\z\ 2k z m 



fc! (fc + to)! 



2k n *m imip\ _ p / \ _ jm(2|jj (Q)) 

g 6 >- Cm ^- /o(2k|<Q)) 



(C8) 



(C9) 



(CIO) 



Expanding this equation to order \z\ 2k z m , one obtains a proportionality relation between the cumulant 

(\Q\ 2k Q* m e lm ^) and (Q) 2k+m v' m . Having measured independently the integrated flow (Q), one thus obtains the 
differential flow v' m from the cumulant. As discussed in Sec. [V, the corresponding error from nonflow correlations is 
of order ((Q) VM)- fe -( m / 2 ). 



4. Removing autocorrelations 



In the case when the "proton" is included in the construction of the event flow vector Q n , i.e. if ip is one of the 
angles (/>,- in Eq. (fl7]), the resulting autocorrelations can be re mov ed at the level of the generating function C m (z) in 
Eq. (C6): this subtraction is similar to that performed in Sec. B4 for the integrated flow. 



Neglecting correlations for simplicity, the generating function of the cumulants, defined by Eq. (|C4|) , becomes 

g(2x cos / ip+2y sin ■ift+imip) / \f~M 



Cm {z) 



g(2x cos ip+2y sin i/>)/%/M \ 



(Cll) 



where we have set z = x + iy, and the brackets denote an average over ip. Assuming for simplicity that the tp 
distribution is isotropic, one obtains 



C m {z) 



I m (2|z|/VM) 
I (2\z\/Vm) \\*\ 



(C12) 



This is the value of the generating function if there are only autocorrelations. If there is flow in addition, we assume 
that the contributions of autocorrelations and flow are additive. Equation (CIO) is then replaced by 



\ Z \ 2k Z m 

^ fc! (fc + m)! 



«|Q| 2fe g* m e m ^)) = I Ir A 2 ^ {Q)) v' 



Io(2\z\(Q)) m h (2\z\/VM 




(C13) 



This equation is equivalent to Eq. (|66|), setting x = \z\. This formula removes exactly all autocorrelations when the 
vector Q n is defined with unit weights. 
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APPENDIX D: INTERPOLATION FORMULAS 



In this Appendix, we give interpolation methods to calculate numerically the cumulants from their generating 
functions. 



1. Integrated flow 

The cumulants used for the measurement of the integrated flow are defined by Eq. (|78|). In order to compute 
numerically the cumulants ((\Q 2k |)) for k' = 1 . . . k from the generating function, one can for instance tabulate the 
generating function at the following points: 

G p ^ q = logt/o ( ^o%/P cos — ^— -,roV"p sin^) (Dl) 

\ ^max (/max / 

for p — 1, . . . , k and q = 0, . . . , <7max — 1. In this equation, ro is a real number which should be chosen small enough for 
the series expansion to converge rapidly, typically ro ~ 0.1, and q max is the number of angles at which the generating 
function is evaluated, which should satisfy the condition g max > 2k. 

One then averages over the angle, thereby eliminating nonisotropic terms up to order \z\ 2k : 

° P = - — (° 2 ) 

4max g=Q 

Then, the G p , with p = 1, . . . , k, are related to the cumulants ((|Q| 2fc )) with k! = 1, . . . , k by the following linear 
system of equations: 

k 2k' 

G P -T,(\Q\ 2k '))^P k ' ^<P<k- (D3) 



For practical purposes, it is enough to take k = 3, as explained in Sec. |IID|. In this case, the solution of the above 
system reads 

{\Q\ 2 }=±(zG 1 - 3 -G 2+ 1 -G 3 

(I0I 4 )) =4 (-5G 1+ 4G 2 - G 3 ), 
r o 

(l<5| 6 » = "I ( 3G i~ 3 G 2 + G 3 ). (D4) 



2. Differential flow 



The cumulants used for the measurement of the harmonic v' m are defined from the generating function by Eq. (|83|). 
In order to compute numerically the cumulants ((|Q| 2fc Q* m e lm ^}} for k' = 0, . . . , k, from the generating function, we 
first tabulate the real and imaginary parts of the generating function, defined by Eq. (183), at the following points: 



Xp, q = 3? 



Y = c > 

1 p.q — ° 



. _ 2qir . 2qw 
C m [ r ^/p cos- , r y/p sin 

QmsiX 

2qir 



?max 

C m I ro^Jp cos^—, r Q ^/p sin 



for p — 1, . . . , k + 1 and q = 0, . . . , 
as we see below. 

One then multiplies C m (z) by z* 



<7max 

— 1. The number of angles q n 



<Zmax y 

x must satisfy the condition q n 



(D5) 
> 2(k + m), 



takes the real part and averages over azimuthal angles. Provided g max is large 



enough, one thus eliminates all nonisotropic terms up to order z' 



k yk-\-m 



in the generating function: 
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a 



m q n 



'In 



E 

9=0 



cos m 



2qir 
I 

?raax 



+ sin ( to 



Then, the values of C p for p = 1, . . . fc + 1 are related to the cumulants 
following linear system of equations: 



2g?r 

9max 
2k' 



Y, 



p.q 



(D6) 

m e""*)) for k' = 0, . . . , k by the 



c P = E((i<3r'Q 



2(fe'+m) fc'+ m 
' n 



fc'!(fc' + to)! 



, 1 <p < fe + 1. 



(D7) 



Taking fc = 1 is sufficient for most purposes, as shown in Sec. IV C. For m = 1, the solution of this system is 



101 



1 



2C 1 - l -C 2 ). 



-L(-2d+C 2 ), 



(D8) 



while for to = 2, 
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4Ci--C 2 ), 
-6d + |c 2 



(D9) 



For a review, see J.-Y. Ollitrault, Nucl. Phys. A638, 195c (1998). 

P. Danielewicz and G. Odyniec, Phys. Lett. 15TB, 146 (1985). 

P. M. Dinh, N. Borghini, and J.-Y. Ollitrault, Phys. Lett. B477, 51 (2000). 

N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, Phys. Rev. C62, 034902 (2000). 

P. Danielewicz et al, Phys. Rev. C38, 120 (1988). 

H. Sorge, Phys. Rev. Lett. 82 (1999) 2048. 

S. A. Voloshin and A. M. Poskanzer, Phys. Lett. B474 (2000) 27. 

L. V. Bravina, A. Faessler, C. Fuchs, and E. E. Zabrodin, Phys. Rev. C61 (2000) 064902. 
A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C58, 1671 (1998). 

M. M. Aggarwal et al, WA93 Collaboration, Phys. Lett. B403 (1997) 390; W.H. van Heeringen, Ph.D. Thesis, ISBN 
90-393-1325-3, Universiteit Utrecht (1996). 



K. H. Ackermann et al. (STAR Collaboration), preprint nucl-ex/0009011 . 

J. Jiang et al, Phys. Rev. Lett. 68, 2739 (1992). 

S. A. Voloshin and Y. Zhang, Z. Phys. C70, 665 (1996). 

J. Barrette et al. (E877 Collaboration), Phys. Rev. Lett. 73, 2532 (1994). 

B. Lenkeit for the CERES Collaboration, Nucl. Phys. A661, 23c (1999). 

A. M. Poskanzer and S. A. Voloshin, unpublished. 

H. Appelshauser et al. (NA49 Collaboration), Phys. Rev. Lett. 80, 4136 (1998). 
J. Barrette et al. (E877 Collaboration), Phys. Rev. C55, 1420 (1997). 
T. L. Hill, Statistical Mechanics, (McGraw-Hill, NY, 1956), chapter 5. 

N. G. van Kampen, Stochastic processes in physics and chemistry (North-Holland, Amsterdam, 1981). 
P. Carruthers and I. Sarcevic, Phys. Rev. Lett. 63, 1562 (1989). 
A. Giovannini and G. Veneziano, Nucl. Phys. B130, 61 (1977). 

H. C. Eggers, P. Lipa, P. Carruthers, and B. Buschbeck, Phys. Lett. B301, 298 (1993). 

J.-Y. Ollitrault, Phys. Rev. D46, 229 (1992); Phys. Rev. D48, 1132 (1993). 

W. K. Wilson, R. Lacey, C. A. Ogilvie, and G. D. Westfall, Phys. Rev. C45, 738 (1992). 

J.-Y. Ollitrault, Nucl. Phys. A590, 561c (1995). 

P. Saturnini (NA50 Collaboration), Nucl. Phys. A661, 345c (1999). 

A. M. Poskanzer and S. A. Voloshin (NA49 Collaboration) Nucl. Phys. A661, 341c (1999). 
P. Danielewicz and M. Gyulassy, Phys. Lett. 129B, 283 (1983). 



34 



[30] H. Gustafsson et al., Phys. Rev. Lett. 52, 1590 (1984). 
[31] K. G. R. Doss et al, Phys. Rev. Lett. 59, 2720 (1987). 

[32] J.-Y. Ollitrault, "A method of reconstructing azimuthal distributions in heavy ion collisions," in XXV International 
Symposium on Multiparticle Dynamics, Stara Lesna, Slovakia, 11-16 September 1995, Bruncko P ., Sandor L., Urban J. 
eds., (World Scientific, 1996), pp. 290-296; J.-Y. Ollitrault, Los Alamos preprint |rmcl-ex/9711003| , unpublished. 

[33] J. Barrette et al. (E877 Collaboration), Phys. Rev. C56, 3254 (1997); Phys. Rev. C59, 884 (1999). 

[34] N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, in preparation. 

[35] C. Domb and M. S. Green Eds., Phase transitions and critical phenomena, vol.3: Series expansions for lattice models 
(Academic Press, New York, 1974). 



35 



